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U.ARSTRACT  (M»V/num  iOQwordi) 

The  simplicity  of  the  pulsed  plasma  micro  thruster  (PPT)  allowed  its  selection  for  space  flight 
application,  even  with  low  efficiency  relative  to  other  concepts.  Improvements  to  PPT  thrust 
efficiency  and  specific  impulse  represent  today’s  challenges  for  research  and  development.  Self- 
consistent  modeling  of  PPT  behavior  requires  a  description  of  the  flow  composition  as  it  changes 
from  complex  molecular  forms  to  highly-ionized,  constituent  atoms.  Heat  flux  from  the  plasma 
discharge  to  the  solid  propellant,  (traditionally,  Teflon),  depends  on  transport  coefficients  for 
electrical  and  thermal  conductivity  in  the  partially  ionized  flow.  Numerical  simulations  now  use 
an  existing,  single-temperature  equation  of  state  for  Teflon  (from  the  SESAME  tables)  along  with 
classical  transport  formulas  based  on  Coulomb  collisions.  The  principal  effort  under  the  present 
grant  has  been  to  develop  a  two-temperature,  JLTE  model  for  Teflon  in  the  regime  of  interest  for 
PPTs.  This  model  includes  25  species  (atoms,  molecules,  ions  and  electrons)  and  allows  separate 
heavy-particle  and  electron  temperatures.  An  idealized  analysis,  limited  to  one-dimensional, 
quasi-steady  MHD  flow,  but  incorporating  resistive  heating  and  thermal  diffusion,  has  also  been 
developed  and  provides  useful  guidance  on  PPT  operation  and  propellant  behavior. 
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INTRODUCTION 


Electric  propulsion  provides  higher  exhaust  speeds  than  chemical  techniques, 
thereby  offering  considerable  economic  advantage  for  several  applications  in  space. 
Today,  these  applications  include  on-board  propulsion  for  satellites  that  already  require 
significant  electrical  power  to  accomplish  their  missions  (e.g.,  communication).  Such 
propulsion  presently  serves  station-keeping  needs,  but  may  also  extend  to  re¬ 
positioning  geosynchronous  satellites,  drag  compensation  and  orbit-raising.  While 
electric  propulsion  no  longer  waits  for  dedicated  power  sources  and  primary  propulsion 
missions,  flight  experience  gained  in  near  earth  missions  may  soon  encourage 
application  of  electric  propulsion  for  planetary  and  deep  space  exploration. 

Basic  research  can  contribute  to  the  application  of  electric  propulsion  by 
providing  the  tools  and  insights  needed  for  improvements  in  specific  impulse,  thrust 
efficiency  and  lifetime.  The  present  research  recognizes  that  transport  processes  play  a 
major  role  in  determining  the  behavior  of  plasma  thrusters.  Modem  computational 
techniques  make  it  possible  to  simulate  complex  plasma  and  electromagnetic 
interactions  within  plasma  thrusters,  but  depend  for  accurate  results  on  the  properties  of 
the  propellant  material.  For  devices  such  as  the  pulsed  plasma  microthruster  (PPT), 
Teflon  has  been  the  propellant  of  choice.  Self-consistent  modeling  of  PPT  behavior 
requires  a  description  of  the  composition  of  Teflon  as  it  changes  from  complex 
molecular  forms  to  highly-ionized,  constituent  atoms.  Knowledge  of  this  composition 
allows  calculation  of  thermodynamic  and  transport  properties.  This  report  describes  the 
development  of  a  model  for  the  composition  of  Teflon  plasma  in  which  separate  electron 
and  heavy-particle  temperatures  exist. 

During  the  period  of  the  grant,  the  Principal  Investigator  had  additional  duties  as 
visiting  Chief  Scientist  for  Advanced  Weapons  and  Survivability,  Phillips  Laboratory, 
Kirtland  AFB,  NM.  Part  of  these  duties  included  examination  of  diffusive  processes  in 
plasma  devices,  such  as  high-power  microwave  sources.  The  particular  problem  of 
impedance  collapse  in  high-voltage,  electron-beam  diodes  involves  plasma  layers  in 
close  contact  with  electrodes.  Modeling  such  layers  with  the  MACH2  code,  a  subject 
related  to  numerical  simulation  of  plasma  thrusters,  led  to  a  technical  paper  included  as 
an  appendix  to  this  report.  Technical  discussions  at  Phillips  Laboratory,  as  part  of  the 
work  on  plasma  closure,  resulted  in  improvements  to  the  treatment  of  both  resistive 
transport  and  non-neutral  plasma  regions  in  the  MACH2  code.  These  improvements 
assist  the  use  of  MACH2  for  modeling  various  kinds  of  electric  thrusters. 


BACKGROUND 

Over  the  past  several  decades,  research  and  development  of  electric  propulsion 
has  embraced  a  great  variety  of  candidate  techniques.  Standard  texts  provide  the 
taxonomy  of  these  techniques.  The  three  main  categories  of  electrothermal, 
electromagnetic  and  electrostatic  propulsion  comprise  dozens  of  concepts.  Plasma 
thrusters  often  combine  aspects  from  more  than  one  category  of  electric  propulsion. 
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They  tend  to  share  characteristics  of  simplicity  and  robustness  of  construction,  and 
typically  have  electrical  impedance  values  in  the  range  of  1  - 100  mQ.  Plasma  thrusters 
operate  in  steady  state  and  pulses  of  duration  down  to  sub-microsecond.  The  inherent 
physics  of  plasma  thruster  operation  may  remain  the  same  over  this  wide  range  of  time 
scales.  If  all  conditions  stay  approximately  constant  for  times  much  longer  than  the  time 
for  flow  transit  through  the  thruster,  operation  is  often  termed  ‘quasi-steady’.  Propellant 
supply,  by  gas  injection  vs  ablation  offers  an  additional  distinction  for  plasma  thrusters. 
Furthermore,  the  magnetic  field  in  the  thruster  may  originate  with  plasma  currents  or  by 
application  of  external  magnets. 

The  present  research  focuses  on  a  form  of  plasma  propulsion  embodied  by  the 
pulsed  plasma  microthruster,  generally  abbreviated  as  PPT.  In  its  traditional, 
rectangular  form,  driven  by  an  LRC-circuit,  a  Teflon  propellant  bar  separates  two 
electrodes  connected  to  a  charged  capacitor  in  vacuum.  A  spark  plug  provides  a  small 
amount  of  initial  plasma,  triggering  an  electrical  discharge  across  the  exposed  surface 
of  the  Teflon.  Heat  transfer  from  this  discharge  causes  evaporation  of  propellant 
material,  which  then  accelerates  through  the  discharge  due  to  electromagnetic  and 
pressure  forces.  As  propellant  evaporates  with  each  discharge,  a  simple  spring 
mechanism  advances  the  propellant  bar  into  the  thrust  chamber.  Two  limiting  modes  of 
operation  can  occur  depending  on  the  details  of  the  heat  transfer  and  acceleration 
process.  For  sufficient  heat  transfer,  the  surface  provides  new  electrically  conducting 
material,  so  the  discharge  path  can  remain  adjacent  to  the  propellant  surface  in  an 
‘ablation  arc’  mode.  Otherwise,  the  discharge  must  follow  the  material  it  accelerates,  so 
the  PPT  operates  in  the  ‘propagating’  mode. 

Choice  of  the  PPT  for  a  research  focus  follows  from  its  special  status  within 
electric  propulsion.  For  many  years,  of  all  the  various  electric  thruster  concepts,  only  the 
PPT  had  achieved  acceptance  for  application  on  actual  space  missions.  Recently,  with 
the  continuing  success  of  the  kilowatt-level  arcjets  on  Telstar  IV,  electric  propulsion 
should  see  representation  as  well  by  xenon  ion  engines  and  Hall  thrusters,  such  as  the 
SPT-100.  In  contrast  with  the  PPT,  these  other  devices  have  all  attained  a  high  degree 
of  refinement  by  decades  of  development  in  the  laboratory,  achieving  perhaps  their 
optimum  levels  of  performance.  For  example,  the  thrust  efficiencies  of  xenon  ion 
engines  and  the  SPT-100  can  exceed  70%,  while  flight  models  of  the  PPT  have 
efficiencies  of  less  than  10%.  The  simplicity  of  the  PPT  allowed  its  selection  for  space 
flight  application,  even  with  such  low  efficiency  relative  to  other  concepts.  Improvements 
to  PPT  thrust  efficiency  and  specific  impulse  represent  today’s  challenges  for  research 
and  development.  Only  the  PPT  combines  decades  of  flight  use  with  considerable 
opportunity  for  increasing  thrust  efficiency  and  specific  impulse. 


APPROACH 

The  experimental  simplicity  of  the  PPT  has  permitted  empirical  studies  to  provide 
sufficient  data  and  experience  to  allow  development  of  space-qualified  hardware  used 
in  actual  missions  without  the  benefit  of  a  complete  theoretical  understanding. 
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Unfortunately,  these  studies  have  failed  to  improve  PPT  efficiency  (with  notable 
exceptions  at  high  energy)  above  several  percent.  Furthermore,  little  insight  has  been 
forthcoming  to  guide  the  choice  of  candidate  propellants  that  might  offer  higher  specific 
impulse  or  better  plume  characteristics  than  the  traditional  use  of  Teflon.  Instead, 
graphs  summarize  empirical  results  for  impulse-bit  and  mass  loss  per  shot  vs  stored 
energy  in  the  PPT  capacitor.  From  these  data,  empirical  constants  permit  closure  of 
simple  models  for  PPT  performance,  thereby  allowing  design  of  ablation-fed,  MPD 
thrusters  .  Minor  manipulation  of  system  equations,  along  with  data  on  masses  and 
lifetimes  of  components  (e.g.,  capacitors),  can  indicate  optimum  parameters  for  PPT 
application  to  particular  missions.  The  renewed  attention  to  the  PPT  for  small  satellites, 
has  created  a  need  for  theoretical  tools  that  can  predict  directions  for  improving  PPT 
performance  in  terms  of  mass  utilization,  efficiency  and  component  reliability. 

The  problem  of  associating  the  mass  evolved  from  the  propellant  surface  with  the 
impulse  developed  during  the  operating  pulse  represents  the  principal  difficulty  in 
establishing  the  performance  of  the  PPT  in  terms  of  its  specific  impulse  or  (average) 
exhaust  speed.  Empirically,  measurement  of  the  total  mass  lost  by  the  surface  over 
many  discharges  provides  the  mass  loss  per  shot.  This  value  divided  into  the  impulse 
per  shot  offers  the  average  exhaust  speed.  Experimental  data  suggest  that  the  mass 
loss  is  proportional  to  the  stored  energy.  This  intuitively  pleasant  result  implies  a 
constant  value  for  the  average  exhaust  speed,  and  a  reasonable  way  to  scale  new 
designs  from  existing  data.  At  this  level  of  discussion,  however,  there  is  no  further 
insight  to  provide  a  firm  basis  for  improving  PPT  performance. 

It  is  necessary  to  develop  theoretical  modeling  that  can  incorporate  physical 
processes  not  captured  by  lumped-circuit  representations  and  empirical  scaling  laws. 
For  example,  the  discharge  current  must  have  the  opportunity  to  flow  both  through 
paths  propagating  along  the  accelerator  electrodes  and  through  a  path  remaining  on  the 
propellant  surface.  This  requires  a  formulation  at  least  at  the  level  of 
magnetohydrodynamics  in  order  to  obtain  the  distribution  of  current  density  within  the 
thruster.  An  effort  to  extend  PPT  operation  to  the  millipound  level  first  attempted  such 
an  approach.  At  the  time,  however,  computational  simulation  of  plasma  acceleration 
was  too  difficult  to  achieve  within  a  program  largely  devoted  to  experimental  exploration 
and  development.  Over  the  past  decade,  modern  calculational  techniques,  created  to 
support  experiments  with  powerful  plasma-guns,  have  been  applied  to  plasma  thruster 
problems.  These  applications  included  quasi-steady  MPD  arcjets  and  steady  state, 
applied-field,  MPD  thrusters,  and  more  recently  PPTs. 

The  need  for  self-consistent  addition  of  mass  during  the  discharge  pulse  makes 
numerical  simulation  much  more  difficult  for  the  PPT  than  fora  gas-fed,  plasma  thruster. 
Presumably,  heat  transfer  from  the  plasma  discharge  to  the  exposed  surface  of  the 
propellant  results  in  evaporation  of  material  through  which  the  discharge  current  flows. 
Electromagnetic  forces  accelerate  this  material,  while  resistive  heating  contributes  to 
the  overall  flow  enthalpy.  Typically,  pulsed  plasma  discharges  propagate  with  the 
accelerated  plasma,  in  order  to  follow  the  electrically  conducting  material,  unless 
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discharge  processes  create  new  conducting  material  closer  to  the  source  of 
electromagnetic  power. 

The  PPT  traditionally  ignites  when  a  spark  plug  generates  an  initial  plasma 
between  the  accelerator  electrodes.  If  new  material  does  not  evaporate  quickly  enough 
from  the  propellant  surface,  the  discharge  will  accelerate  this  initial  plasma,  travelling 
with  it  along  the  electrodes.  The  increased  separation  of  the  discharge  from  the 
propellant  surface  reduces  the  opportunity  for  heat  transfer  and  thus  encourages  a 
propagating  mode  for  discharge  operation.  Similar  behavior  would  occur  if  the  discharge 
current  waveform  resulted  in  rapid  acceleration  of  an  initial  amount  of  material 
evaporated  from  the  propellant  surface  before  heat  transfer  was  sufficient  to  maintain 
further  mass  addition.  On  the  other  hand,  adequate  heat  transfer  to  the  surface  can 
provide  a  continual  source  of  electrically  conducting  material,  allowing  the  discharge  to 
remain  adjacent  to  the  surface,  thereby  sustaining  heat  transfer. 

To  determine  the  actual  mode  of  operation  of  the  PPT  for  various  possible  values 
of  parameters  and  arrangements  requires  time-dependent  calculations  combining 
electrical  circuitry,  MHD  flow  and  heat  transfer  in  the  propellant.  The  MACH2  code, 
modified  to  include  the  features  needed  for  simulating  the  PPT,  was  first  used  to  model 
the  LES-6  thruster.  Modifications  included  the  development  of  a  separate  numerical 
model  for  two-dimensional,  unsteady  heat-flow  in  the  solid  propellant,  based  on  the  net 
heat  flux  to  the  exposed  surface.  This  heat  flux  comprised  heat  conduction  and  radiation 
from  the  plasma  discharge,  and  convection  due  to  evaporation  or  condensation  of  the 
propellant  material.  (Earlier  models  for  ablation  of  solids  assumed  a  thermal  diffusion 
depth  that  increased  with  the  square  root  of  elapsed  time.  Such  an  approach  might 
succeed  for  monotonically  increasing  temperatures,  but  cannot  service  the  variations 
possible  in  the  PPT.)  The  initial  calculations  for  LES-6  captured  both  the  magnitude  and 
variation  of  the  impulse-bit  for  the  available  values  of  experimental  data.  Such 
agreement  must  be  considered  somewhat  fortuitous,  the  result  of  compensating  errors, 
because  of  several  limitations.  These  include  the  use  of  a  two-dimensional  calculation 
(in  the  plane  of  the  current  flow)  for  a  three-dimensional  problem,  and  approximate 
formulations  for  plasma  composition  and  transport  coefficients. 

Heat  flux  from  the  plasma  discharge  depends  on  the  distribution  of  current 
density  near  the  propellant  surface  and  the  transport  coefficients  for  electrical  and 
thermal  conductivity  in  the  partially  ionized  flow.  Present  calculations  use  an  existing, 
single-temperature  equation  of  state  for  Teflon  (from  the  SESAME  tables)  along  with 
classical  transport  formulas  based  on  Coulomb  collisions.  The  principal  effort  under  the 
present  grant  has  been  to  develop  a  two-temperature,  LTE  model  for  Teflon  in  the 
regime  of  interest  for  PPTs.  This  model  includes  25  species  (atoms,  molecules,  ions 
and  electrons)  and  allows  separate  heavy-particle  and  electron  temperatures. 

An  idealized  analysis,  limited  to  one-dimensional,  quasi-steady  MHD  flow,  but 
incorporating  resistive  heating  and  thermal  diffusion,  has  also  been  developed  to 
provide  some  guidance  on  PPT  operation  and  propellant  behavior.  In  particular, 
application  of  a  magneto-sonic  choking  condition  at  the  downstream  edge  of  the 
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discharge  determines  the  mass  flow  rate  needed  by  the  MHD  flow  for  given  values  of 
total  current  and  size.  This  mass  flow  rate  in  turn  specifies  the  surface  temperature  of 
the  propellant  and  the  associated  equilibrium  vapor  pressure.  The  heat  delivered  to  the 
propellant  to  achieve  these  conditions  then  depends  on  the  given  current  waveform  and 
duration. 


PROGRESS 

The  five  quarters  of  sponsorship  under  the  present  grant  saw  accomplishment  of 
the  first  portions  of  the  proposed  two-year  effort.  Appendix  I  displays  much  of  this  work 
in  the  form  of  a  presentation  given  at  the  AFOSR  workshop  in  San  Diego,  28  -  31  July 
1997.  Appendices  II  and  III  provide  copies  of  technical  papers  on  the  directions  for  PPT 
improvement  (based  on  the  idealized  model)  and  the  two-temperature,  LTE  model  for 
Teflon.  The  25th  International  Electric  Propulsion  Conference,  Cleveland,  OH,  24  -  28 
August  1997,  served  as  the  forum  for  these  papers,  which  will  appear  in  the  conference 
proceedings.  Appendix  IV  provides  an  updated  version  of  the  paper  on  transport 
properties  of  nitrogen,  which  led  to  the  concern  with  severely  non-monotonic  variations 
of  thermal  conductivity  with  temperature  and  pressure  in  molecular  gases.  This  concern 
resulted  in  attention  (Appendix  V)  to  capturing  such  variations  in  numerical  simulations 
by  means  of  computational  grids  that  adapt  to  thermal  conductivity,  rather  than  simply 
geometry  or  flow  density.  Appendix  VI  presents  the  latest  form  of  the  two-temperature, 
LTE  model  for  Teflon,  including  more  accurate  calculation  of  the  effects  of  rotation  and 
vibration  of  polyatomic  species.  The  new  model  also  allows  for  variation  of  coupling  of 
vibrational  states  to  electrons  vs  heavy-particles. 

Appendix  VII  provides  a  paper,  presented  at  the  11th  IEEE  International  Pulsed 
Power  Conference,  Baltimore,  MD,  30  June  -  3  July  1997,  on  impedance  collapse  in 
high-voltage,  electron-beam  diodes.  It  summarizes  the  work  on  diffusive  processes  near 
electrodes,  performed  by  the  Principal  Investigator  while  at  Phillips  Laboratory,  Kirtland 
AFB. 


CONCLUDING  REMARKS 

The  most  recent  theoretical  activity  at  Ohio  State  in  support  of  PPT  development 
has  been  under  NASA  sponsorship  and  comprises  continuation  of  MACH2  calculations 
applied  to  the  benchmark  PPT,  and  extension  of  such  calculations  to  coaxial 
configurations.  The  latter  effort  draws  on  earlier  success  in  using  MACH2  to  model 
quasi-steady  MPD  thrusters.  Simulations  of  coaxial  PPTs  have  already  included  the  use 
of  plug  nozzles  to  provide  improved  flow  expansion.  Modification  of  MACH2  to  couple 
the  PPT  to  other  circuits,  such  as  PFNs  and  the  inductively-driven  circuit  has  also 
begun. 


Under  AFOSR  sponsorship,  transport  coefficients  based  on  the  two-temperature 
LTE  model  for  Teflon  are  now  developing.  Extension  of  such  modeling  to  other 
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candidate  materials  could  guide  propellant  selection.  Decomposition  of  solid  propellants 
exposed  to  high  current  discharges  remains  a  critical  area  of  concern  for  both  propellant 
choice  and  theoretical  modeling.  The  idealized  model  has  already  suggested  a 
connection  between  discharge  operation  and  material  properties.  At  sufficiently  high 
currents,  the  temperature  needed  to  supply  the  mass  flow  to  the  discharge  will  exceed 
values  at  which  the  solid  propellant  decomposes  in  some  manner.  For  example,  the 
surface  may  liquefy.  A  thin  liquid  film  pressed  by  vapor  against  the  still  solid  interior  of 
the  propellant  can  splash  laterally  on  the  electrode  and  insulator  surfaces  of  the  thrust 
chamber,  resulting  in  mass  ejection  at  low  speed.  If  such  a  mechanism  exists  in  present 
PPTs,  it  can  provide  a  significant  limitation  on  PPT  performance.  Proper  choice  of 
propellant  and  operating  magnetic  field  to  avoid  pulsed  liquid  films  would  avoid  this 
difficulty.  The  results  of  the  idealized  model  depend  on  transport  properties,  such  as  the 
thermal  conductivity.  Decomposition  of  solid  propellants  involves  particular  values  of 
temperature  (e.g.,  melting  point).  Combination  of  plasma  transport  and  solid 
decomposition  will  require  self-consistent  calculations  based  on  accurate  plasma 
properties  and  careful  numerical  modeling. 
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ABSTRACT 

Potential  improvements  in  PPT  performance 
are  discussed  from  analytical  considerations,  and 
comprise  the  following  sequence  of  conditions  that 
need  to  be  established.  Presently,  the  exhaust  of  the 
PPT  is  allowed  to  expand  without  regard  to  extracting 
directed  kinetic  energy  efficiently  from  die  hot, 
highly-magnetized  plasma  in  the  thrust  chamber. 
Analytically,  it  appears  that  improvements  by  a  factor 
of  1.7  in  specific  impulse  and  three  in  thrust 
efficiency  should  be  possible  with  proper  expansion. 
A  major  area  for  improving  PPT  performance  is  the 
reduction  of  relative  mass  expelled  that  is  not 
electromagnetically-accelerated  to  high  speed. 
Analytical  modeling  indicates  that  mass  evolved  by 
post-discharge  evaporation  can  exceed  that  during  the 
discharge  by  a  factor  of  more  than  five.  An 
inductively-driven  circuit,  previously  suggested, 
would  maintain  electromagnetic  acceleration  as  the 
propellant  surface  cools.  This  circuit  also  improves 
PPT  performance  by  eliminating  losses  and 
difficulties  of  present  oscillatory  waveforms. 

INTRODUCTION 

It  is  useful  to  consider  directions  for 
improving  the  pulsed  plasma  microthruster  (PPT)  so 
that  it  may  be  applied  to  a  greater  range  of  missions. 
In  particular,  higher  thrust  efficiency,  and,  for  some 
uses,  higher  specific  impulse  are  needed.  It  is  critical 
that  improvements  to  the  PPT  retain  the  simplicity 
that  allowed  its  early  operation  in  actual  space 
missions,  and  maintain  the  connection  to  the  PPTs 
extensive  flight-experience.  Of  the  four  or  five 
devices  selected  from  dozens  of  concepts  in  electric 
propulsion,  only  the  PPT  has  both  a  record  of  actual 
accomplishment  in  space,  mid  the  potential  for 
significant  improvement  through  further  research. 
The  other  devices  have  already  been  taken  to  high 
levels  of  performance  (probably  their  limiting  values) 
by  many  years  of  sustained  laboratory  research  and 
development 
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The  traditional  pulsed  plasma  microthruster 
(PPT)  operates  with  an  unsteady,  oscillatory 
discharge.  This  is  a  consequence  of  the  relatively  low 
energies  used  by  PPTs  for  satellite  station-keeping. 
With  a  stored-energy  in  the  capacitor  on  the  order  of 
20  J,  at  an  initial  voltage  of  2  kV,  the  capacitance  is 
only  10  pf.;  circuit  inductance  values  below  100  nh 
are  difficult  to  obtain  with  commercially-available 
components.  The  impedance  of  an  LRC -circuit  for 

critical  damping  is  2(L/C)^  =  200  mfl  The 
characteristic  impedance  of  the  electromagnetically- 
accelerated  discharge  flow,  however,  is  (for  a 
propagating  discharge,  with  inductance  gradient,  L/): 

Zd  =  L'u/2  (1) 

so,  at  an  exhaust  speed  of  40  km/s,  and  L'  = 
10~6  h/m,  Zd  =  20  m£l.  Thus,  the  circuit  is  hardly 
loaded  by  the  thruster,  and  can  deposit  much  of  its 
energy  in  the  internal  resistance  of  the  capacitor.  The 
losses  associated  with  such  resistance  increases  the 
operating  temperature  of  the  capcacitor.  The 
combination  of  higher  temperature  with  the  severe 
and  repeated  reversals  of  the  capacitor  voltage 
reduces  the  reliable  life  of  the  capacitor,  which  must 
be  compensated  by  a  reduction  in  design  voltage  and 
energy  per  unit  mass. 

It  is  also  typical  of  traditional  PPT  operation 
that  the  discharge  flow  simply  exits  abruptly  from  a 
constant  area.  For  a  magnetized-plasma  flow,  this 
fails  to  extract  energy  from  the  magnetic  field  (and 
any  accessible  thermal  modes)  into  the  directed 
kinetic  energy  of  the  exhaust.  An  additional  source  of 
inefficiency  in  traditional  PPTs  has  recently  been 
suggested  after  initial  attempts  at  numerical 
simulation  of  the  LES-6  device  The  numerical 
calculations,  using  the  MACH2  code,  agreed  well 
with  the  impulse-bit  per  shot,  but  the  mass  ablated 
during  the  discharge  pulse  was  about  a  factor  of  ten 
lower  in  the  simulation  vs  the  experimental  data  for 
mass  loss  per  shot.  Subsequent  calculation  indicated 
that  evaporation  of  the  Teflon  propellant  between 
firings  might  account  for  this  mass  discrepancy.  If  a 
major  portion  of  the  mass  per  shot  is  lost  at  relatively 
low  speed,  then  the  efficiency  of  the  PPT  is 
substantially  reduced  from  ideal  values. 

In  combination,  die  efficiency  factors 
associated  with  external-circuit  resistance  (<60% ), 
improper  flow  expansion  (<33%)  and  mass  loss  at 
low  speed  (<50%)  multiply  to  provide  a  total  thruster 
efficiency  of  less  than  10%  .  By  addressing  each  of 
the  inefficiencies  in  turn,  it  should  be  possible  to 
improve  the  performance  of  the  PPT  substantially. 


Consideration  may  be  framed  first  in  terms  of  quasi- 
analytical  modeling,  before  invoking  more  powerful 
numerical  tools  and  testing  ideas  experimentally. 
Idealized  modeling  may  be  applied  to  initial  design  of 
proper  flow  expansion,  to  the  questions  of  mass  loss 
and  thermal  management  and,  finally,  to  improved 
circuitry  for  the  PPT. 

IDEALIZED  MODEL  FOR  PPT 

The  essential  features  of  the  ablation-fed 
discharge  in  the  PPT  include  resistive  heating  near 
the  entrance  of  a  constant-area  channel  (where  the 
back  EMF  is  relatively  low),  heat  transfer  from  the 
discharge  back  to  the  propellant  surface  to  provide 
mass  by  ablation,  and  electromagnetic  acceleration  of 
the  plasma  by  the  Lorentz  force.  In  pulsed  operation, 
particularly  with  high-frequency,  oscillatory 
waveforms,  the  preceding  features  are  unsteady  and 
require  numerical  modeling  for  accurate  calculations 
in  time  and  space.  A  first  step  in  simplifying  analysis 
of  the  PPT,  while  attempting  to  retain  the 
fundamental  interactions  among  resistive  heating, 
heat  transfer  and  flow,  restricts  examination  to  a 
steady  state,  and  one-dimensional  flow.  (The  details 
of  such  examination  are  briefly  described  in 
Appendix  L)  The  use  of  a  steady  analysis  in 
discussing  the  PPT,  however,  means  that  comments 
Can  only  be  applied  to  situations  in  which  there  is 
enough  time  for  the  discharge  flow  to  operate  with  a 
balance  of  heat  conduction,  resistive  dissipation  and 
flow  acceleration.  Convective  times  based  on  the 
discharge  thickness  divided  by  the  flow  speed  must 
certainly  be  less  than  the  time  for  variation  of  circuit 
current.  Furthermore,  the  ablating  surface  must  be 
able  to  supply  new  material  in  times  shorter  than  the 
convective  time.  Thus,  for  example,  discharges  that 
lift  off  of  refractory  insulators  may  remain  in  an 
unsteady,  propagating  mode,  rather  than  achieving 
the  quasi-steady  situation  of  the  present  analysis. 

For  a  Teflon-based  PPT,  the  analysis  of 
Appendix  I  suggests  that  an  ablation  arc,  with  a 
thickness  less  than  two  millimeters,  is  formed 
adjacent  to  the  propellant  surface.  In  the  numerical 
example,  the  speed  of  the  flow  through  the  arc 
increases  by  a  factor  of  about  three  to  an  exit  velocity 
of  41  km/s.  The  characteristic  convective  time  for  the 
flow  structure  is  therefore  about  0.1  psec.  Heat 
conduction  to  the  colder  upstream  boundary 
automatically  supplies  the  power  needed  to  dissociate 
and  ionize  the  flow,  and  also  the  relatively  minor, 
additional  power  level  required  to  provide  mass  flow 
by  ablation.  The  calculated  size  and  timescale 
suggest  that  the  analysis  is  consistent  with  PPT 
operation  at  frequencies  (within  the  discharge  pulse) 
less  than  a  MHz,  and  dimensions  greater  than  a  cm. 


SPECIFIC  IMPULSE 

As  previously  noted  in  a  simpler  analysis 
(without  heat  conduction)^,  the  proportionality  of 
resistive  heating  and  electromagnetic  work,  in  the 
context  of  a  flow  in  which  heat  is  largely  absorbed  by 
the  ionization  of  the  propellant,  leads  to  an  exhaust 
speed  that  scales  closely  with  Alfven  critical  speed. 
Thus,  the  specific  impulse  of  self-field,  plasma 
thrusters,  operating  with  mass  addition  (vs  constant 
mass,  propagating  discharges)  will  tend  to  values 
proportional  to  Aliven  critical  speed,  if  heating  can 
supply  additional  conducting  material. 
Improvements  of  PPT  performance,  in  terms  of 
higher  specific  impulse,  therefore,  would  require 
propellants  with  lower  average  molecular-mass  than 
the  Teflon  presently  used.  In  the  present  example, 
the  computed  exit  speed  already  corresponds  to  a 
specific  impulse  of  4180  s.  Even  moderate  attention 
to  proper  expansion  of  this  magnetosonic  flow  to 
magnetic  field-free  conditions  will  offer  values  of 
specific  impulse  that  cover  die  range  of  any  near  term 
missions,  (upwards  of  7000  s).  The  earlier  analysis^ 
suggests  improvements  by  up  to  3^.  The  principal 
reason  for  the  more  modest  values  of  specific  impulse 
(~  1000  s)  is  mass  that  is  not  accelerated  electro- 
magnetically  (e.g.,  post-discharge  evaporation). 

MASS  ABLATED  DURING  DISCHARGE 

By  the  analysis  of  Appendix  I,  the  necessary 
mass-flow  rate  is  actually  controlled  by  a 
magnetosonic  condition  in  the  constant-area  channel, 
rather  than  a  separate  condition  on  heat  transfer  to  the 
propellant  surface.  The  details  of  such  heat  transfer 
adjust  to  satisfy  the  mass  flow  constraints  in  the 
overall  MHD  flow.  The  mass  loss  during  the 
discharge  pulse  may  thus  be  estimated  from  the  mass 
flow  rate  (per  unit  area): 

*  * 

w  =  p  u  (2) 

where  the  speed  at  the  sonic  point  is: 

u*  =  {(B*2  /  p*M)[  1  +  Y0*  /  2  ]} 1/2  (3) 

and  the  mass  density  there  is: 

p*  =  (p*B*2  /  2jj)  /  R*T*  (4) 

The  temperature  at  the  sonic  point  is  obtained  in 
terms  of  the  magnetic  field  at  the  propellant  surface, 
B|: 

T*  =  {A/fi  pflth/Kr)172}2/5  Bi2^  (5) 


The  mass  flow  rate  is  therefore  proportional  to 
For  a  constant  current,  the  mass  ablated  during  a 
pulsetime  tp  is  merely  wAtp,  where  A  is  the  area  of 
the  ablating  surface.  In  the  case  of  an  exponentially- 
decaying  sinusoidal  pulse,  within  the  quasi-steady 
approximation,  the  mass  flow  rate  may  be  integrated 
over  the  oscillatory  waveform;  if  the  ratio  of  risetime 
to  decay  time  is  0.3,  for  example,  the  mass  ablated  is 
0.933woAtr  ,  where  tr  is  the  risetime  and  wG  is  the 
mass  flow  rate  per  unit  area  based  on  the  undamped 
amplitude  of  the  current.  Note  that  this  represents  a 
nearly  linear  dependence  on  stored  energy,  WG,  in  the 
capacitor  (Am  ~  W0^l°)  for  the  mass  ablated  during 
the  discharge. 

THERMAL  CONDITKXNS  AT  SURFACE 

The  mass  lost  between  discharges  may  be 
considered  in  terms  of  the  temperature  of  the 
propellant  surface.  From  the  idealized  analysis,  it  is 
possible  to  estimate  the  surface  temperature  of  the 
propellant  that  is  consistent  with  the  flow  conditions. 
In  particular,  the  equibrium  vapor  pressure  should 
equal  the  total  pressure  at  the  entrance  to  the  ablation 
arc.  A  formula  for  the  equilibrium  vapor  pressure  for 
Teflon  is: 

Peq  =  Pc  exP  (-  Tc  /  T)  (6) 

with  Tc  =  20,815  K  and  pc  =  1.872  x  1020 

N/m2. 

The  total  pressure  calculated  from  the  one¬ 
dimensional,  idealized  model  is: 

pt=pi  +  piui2  (?) 

—  (B i2 / 2jj fj2 HP*  0i/(Oi 

+  (2  +  yP*)031  ] 

The  necessary  surface  temperature  is  then: 

Ts  =  Tc  /  In  {  pc  /  (B  i2  /  2|i  fi2 )[  f  0i/«  l 

+  (2  +  yp*)tfl]}  (8) 

For  the  numerical  example  of  Appendix  I,  the  surface 
temperature  is  600  K.  Note:  this  value  is  only  weakly 
dependent  on  the  operating  magnetic  field,  but 
material  transitions  can  be  quite  sensitive  to  exact 
values  of  temperature.  This  particular  value  is  very 
close  to  the  melting  point  of  Teflon  (~  600  K). 
(Nonuniformities  in  arc  distribution  across  the  face  of 
the  propellant  might  cause  local  melting  in  any 
event.)  The  depth  of  propellant  heated  to  this 
temperature  during  pulsetimes  of  several 
microseconds  is  less  than  a  few  microns.  Growth  of 


perturbations  of  a  liquefied  surface  due  to  Rayleigh- 
Taylor  instability  would  be  suppressed  for 
wavelengths  approaching  the  depth  of  the  layer, 
while  the  exponential  growth  of  shorter  wavelengths 
would  not  persist  beyond  amplitudes  comparable  to 
these  wavelengths.  Thus,  micron-size  droplets  might 
be  expected,  especially  from  edges.  Such  droplets 
would  be  responsible  for  mass  loss  by  surface 
disruption,  as  indicated  in  some  experiments^. 

THERMAL  MANAGEMENT 

It  has  been  suggested^  that  the  loss  of  mass 
between  shots  depends  critically  on  the  overall 
thermal  management  of  the  PPT,  both  in  the 
laboratory  and  in  space.  The  estimated  surface 
temperature  is  well  above  mean-values  within  the 
propellant  measured  in  laboratory  tests^  at  total 
power  levels  of  40  W  (40  J  at  1  Hz)  which  indicate  a 
rise  over  several  thousand  shots  from  room 
temperature  (300  K)  to  about  370  K.  For  the 
acknowledged  low  efficiency  of  present  PPT 
operation,  only  a  small  fraction  of  the  total  power  is 
delivered  to  the  Teflon  surface.  The  estimated  surface 
temperature  allows  calculation  of  the  heat  deposited 
in  the  surface  during  the  pulse,  based  on  the  thermal 
skin-depth ,  6: 

H  =  pc  A5  ( Ts  -  Tj )  (9) 

where  the  thermal  skin-depth  is  given  in 
terms  of  the  pulsetime  tp  as: 

5  =  ( k  tp  /  pc  )l/2  (10) 

With  k  =  0.305  W/m-s-K,  p  =  2. 15  x  103  kg/m3,  and 
c  ss  1171  J/kg,  a  pulsetime  of  10  jjs  would  provide  a 
skin-depth  of  1.1  microns.  At  a  surface  temperature 
of  600  K,  this  contains  637  J/m^  of  heat  added  by  the 
discharge  pulse,  which  represents  an  average  heat 
load  to  a  4  cm^  surface  of  0.25  W  at  a  1  Hz  repetition 
rate.  This  exceeds  the  power  required  to  evaporate 
Teflon  from  the  surface  by  a  factor  of  5  (using  the 
mass  flow  computed  in  Appendix  I  and  a  heat  of 
vaporization  and  de-polymerization  of  3.67  MJ/kg). 
The  "extra"  power  delivered  to  the  surface  has 
consequences  for  both  late-time  pulsed  and  steady 
mass  evolution. 

After  the  discharge  pulse  ends,  the  heat 
deposited  in  the  skin-layer  will  be  shared  with  the  rest 
of  the  solid  propellant  in  a  depth  that  continues  to 
increase  as  the  square-root  of  time.  Without  further 
heat  addition  (or  significant  cooling  due  to  ablation), 
the  surface  temperature  will  decrease  inversely  with 
this  depth: 

(Ts-Tb)/(Tsi-Tb)  =  (tp/t)l/2  (11) 


where  Tsj  is  the  surface  temperature  at  the 
end  of  the  discharge  pulse  (t  =  tp),  and  Tb  is  the  base 
temperature  of  the  propellant.  The  mass  evaporated 
as  the  surface  cools  (for  t  >  tp  )  may  be  estimated 
using  this  time-dependence  of  the  surface 
temperature  in  Eqn.  6. 

By  assuming  one-dimensional  expansion  of 
the  surface  vapor  to  a  (thermally)  sonic  condition,  the 
mass  flow  per  unit  area  is: 

w  =  (Y  /  2)[2/(y  +  1)  ] 1/2  ps  /  ( y  RTS ) 1/2  (12) 

Integration  of  this  mass  flow  rate  provides  a  total 
evaporated  mass  (for  t  >  tp),  Arne,  that  is  proportional 
to  the  magnetic  pressure  and  the  pulsetime,  allowing 
comparison  with  the  mass  ablated.  Any  ,  during  the 
discharge  pulsetime: 

Ame/  Any  >  K{  1  -  exp  (-  Tc  /  Tb)  /  exp  (-  Tc  /  Tsj)} 

(Tc/T*)1/2(TC/Tsi)1/2  (1- 
Tb  /  Tsi)[p*(l  +  Yp*/2)(y  +  1)/  4y  ]  ^ 

(B) 

For  the  conditions  of  the  previous  numerical 
example,  a  base  temperature  Tb  ~  370  K,  and  the 
factor  K  =  1.26,  the  ratio  of  mass  evaporated  as  the 
surface  cools  to  that  ablated  during  the  discharge 
pulse  is  5.1  .  This  ratio  decreases  to  4.3,  if  the  base 
temperature  of  the  propellant  is  kept  at  300  K, 
indicating  improved  performance  of  PPTs  with  better 
cooling.  The  ratio  increases  largely  as  E\  due  to 

the  variation  of  T*,  and  thus  is  rather  insensitive  to 
the  amplitude  of  the  circuit  current.  Longer 
pulsetimes  increase  the  mass  ablated  during  the 
discharge,  but  also  increase  the  heat  deposited  in  the 
solid  propellant,  which  maintains  the  surface 
temperature  for  a  longer  time  after  the  pulse,  allowing 
significant  evaporation  to  continue  longer.  Major 
improvements  will  require  either  optimization  of 
material  properties  or  matching  of  the  power  circuit 
to  the  ablation  process  to  avoid  evolution  of  mass 
when  electromagnetic  forces  are  absent. 

QUASI-STEADY,  INDUCTIVE  OPERATION 

One  approach  to  preventing  such  evolution 
would  reduce  the  surface  temperature  needed  to 
support  the  mass  flow  through  the  discharge  as 
thermal  conduction  into  the  solid  cools  the  surface. 
After  an  initial  pulsetime  due  to  the  current  rise,  tp  = 
tr,  the  surface  temperature  would  then  decline 
according  to  Eqn.  11,  if  no  significant  additional  heat 
deposition  is  required  in  order  to  supply  mass  flow  to 
the  discharge.  Now,  let  the  current  decrease  so  that 


the  required  temperature  follows  the  decreasing 
temperature  of  the  surface: 

J/Jo  =  Bi/Bio 

=  exp  [-  (Tc  /  Tb)  /  (1  +  (Tsi  /  Tb  -  l)(tp  / 1  )1/2] 
/exp  [-  (Tc  / Tb)  / (1  +  (Tsi  / Tb  - 1)] 

(14) 

Hus  is  displayed  in  Fig.  1  with  t  in  units  of  the 
risetime,  tj-  =  tp.  (Temperature  values  are  the  same  as 
in  the  earlier  discussion.)  Such  a  waveform  may  be 
compared  with  the  experimental  current  behavior 
(Fig  2)  obtained  with  an  inductively-driven  circuit^  in 
which  a  plasma  discharge  (in  this  case  a  second  PPT) 
is  used  to  crowbar  the  capacitor  shortly  after  peak 
current. 


Figure  1:  Normalized  current  waveform,  J  /  Jo  vs 
time  in  units  of  risetime,  t  /  tr 


lusec 

Figure  2:  Experimental  current  waveform  for 
inductively-driven  circuit  driving  PPTs 


The  use  of  a  inductively-driven  circuit  not 
only  provides  a  current  waveform  that  might  alleviate 
mass  evolution  after  the  current  pulse,  but  also  avoids 
the  severe  voltage-reversals  on  die  capacitor.  With 
some  attention  to  reducing  the  resistance  of  the 
external  circuit  to  a  small  fraction  of  the  PPT 
impedance,  the  electrical  efficiency  should  greatly 
improve.  From  the  idealized  analysis,  the  impedance 
of  the  PPT  is: 

Z  =  u*B*h/J  (15) 

where  h  is  the  length  of  the  discharge.  For 
the  values  previously  used,  and  h  =  2  cm,  Z  =  33  mQ. 
At  J0  =  10  kA,  and  an  initial  circuit  energy  of  20  J, 
the  inductance  of  the  store  could  be  400  nh,  for  which 
the  characteristic  decay  time  of  the  waveform  would 
be  12  jisec;  the  capacitance  at  an  initial  voltage  of  2 
kV  is  10  jif,  so  the  risetime  is  about  3.1  jisec. 

CONCLUDING  REMARKS 

The  idealized  analysis  has  indicated  that 
evolution  of  mass  after  the  discharge  pulse  is  a 
fundamental  consequence  of  creating  mass  by 
ablation  during  the  discharge.  It  is  therefore  useful  to 
maintain  electromagnetic  forces  while  the  surface 
cools.  This  can  be  accomplished  simply  by  means  of 
an  inductively-driven  circuit,  which  merely  involves 
placing  a  low  impedance  across  die  energy-storage 
capacitor  shortly  after  peak  current  Such  a  circuit 
was  originally  suggested^  to  improve  PPT  design  by 
allowing  high  energy  per  unit  mass  at  low  total 
energies  (without  the  difficulties  of  parasitic 
inductance  and  internal  resistance  in  the  capacitor).  In 
addition  to  reducing  internal  losses,  reduction  of  the 
amplitude  of  voltage-reversal  on  the  capacitor 
improves  reliability  at  high  energy  density. 
Furthermore,  die  new  circuit  provides  longer 
discharge  times,  so  that  proper  flow  expansion 
techniques  can  be  used;  nozzle  sizes  divided  by  flow 
speeds  require  quasi-steady  currents  for  several  jisec. 

While  the  idealized  analysis  can  guide 
general  considerations,  and  may  closely  match 
experimental  data  in  some  cases,  accurate  analysis 
requires  numerical  tools,  such  as  MACH2.  This 
includes  design  of  a  properly  expanded  PPT  flow, 
which  has  recently  been  successfully  attempted  with 
an  annular  PPT  exiting  to  a  plug  nozzle. 
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APPENDIX  I 

ONE-DIMENSIONAL,  STEADY 
PLASMA  THRUSTER  FLOW  WITH 
HEAT  CONDUCTION 

The  conservation  equations  for  a  one- 
dimensional,  steady  plasma  flow  with  resistive 
heating  and  heat  conduction  are  written  as: 

Mass-flow : 

pu  =  constant  =  w  (A  1) 

where  p  is  the  local  mass  density  and  u  is  the 
local  flow  speed.  The  constant  may  be  evaluated  in 
terms  of  values  at  a  particular  location,  w  =  p  u  . 

Momentum  : 

pu  du  +  d_(  B2/2jj  +  p )  =  0  (A2) 

dx  dx 

so,  wu  +  B2/2ji  +p  =  constant 

=  p  V  2  +  B*2/2p  +  p*  (A3) 

where  B  is  the  local  magnetic  field  and  p  is 
the  local  pressure;  starred  quantities  are  evaluated  at 
the  same  location.  For  simplicity,  the  pressure  may  be 
written  as: 

p  =  pRT  (A4) 

where  T  is  the  local  temperature,  and  R  is  an 
appropriate  gas  constant.  The  momentum  equation 
then  becomes: 


wu  +  B2/2p  +  wRT/u 

=  pV  2  +  B*2/2p  +  p*R*T*  (A5) 

The  temperature  distribution  depends  on  the 
interplay  of  convection,  heat  conduction,  resistive 
dissipation,  and  work,  which  may  be  written  in  terms 
of  the  third  conservation  equation: 

Energy : 

w  dU  =  d  k  dT  +  qj^  -  p  du  (A6) 

dx  dx  dx  dx 

where  U  is  the  energy  per  unit  mass,  k  is  the 
thermal  conductivity,  q  is  the  electrical  resistivity, 
and  j  is  the  current  density.  In  steady  state,  and  one- 
dimension,  the  resistive  dissipation  may  be  written  in 
terms  of  an  electric  field  that  is  uniform: 

E  =  qj  -uxB  +(jxB-gradpb)/nee  (A7) 

=  constant  =  E*  (A8) 

where  the  electric  field  can  be  evaluated  at 
the  starred  location  for  which  j  =  j*  ,  so  E*  =  q*j*  + 
u*B*.  (Note  that  the  gradient  of  the  electron 
pressure,  pe,  divided  by  the  electron  density,  does 

not  contribute  to  qj2  in  the  one-dimensional  problem, 
nor  does  the  Hall  effect  term.)  The  energy  equation  is 
then  given  by: 

wdU  =  dkdT  +  (E*-  uB)  -  pdu  (A9) 

dx  dx  dx  q  dx 

In  general,  solution  of  this  equation  can  be 
accomplished  if  the  detailed  behaviors  of  the  energy 
and  pressure  functions,  and  the  transport  coefficients, 
k  and  q,  are  known  in  terms  of  temperature  and 
density. 

For  ideal  plasmas  A 1,  the  thermal 
conductivity  and  electrical  resistivity  can  be  written 
in  terms  of  formulas  that  only  involve  the 
temperature,  if  the  plasma  is  sufficiently  ionized  (and 
does  not  change  its  degree  of  ionization)  and 
magnetic  fields  do  not  suppress  the  electron  heat 
conduction  unduly.  Thus, 

k  =  KhT5/2  and  p  =  Kr/T3/2  (A10) 

where  Kfr  and  K*  are  constants.  (We  have 
also  ignored  the  variation  of  the  Coulomb  logarithm 
here  for  continued  simplicity.) 

To  delineate  the  flow  structure,  and  avoid 
losing  general  results  in  consideration  of  particular 
plasma  values,  it  is  useful  to  non-dimensionalize 
variables  in  terms  of  conditions  at  the  starred  location 


(which  might  later  be  identified  as  a  sonic  point). 
Thus,  let: 

6  =  T/T*,  eo  =  u  /  u*  •  f  =  B  /  B*  , 

and  a  =  x / Xc  (All) 

where  Xc  is  a  characteristic  distance 
determined  later.  Three  dimensionless  parameters  are 
also  obtained: 

=  p  / (B  - / 2p  ) ,  Rm  =  u  B  /p  j 

and  P  =  wcpxc  /  KhT*5/2  (A  1 2) 

where  p*  and  Rm*  are  the  plasma-beta  and 
local  magnetic  Reynolds  number  at  the  starred 
location,  and  P  is  essentially  a  Peclet  number  based 
on  the  characteristic  length: 

Xc  =  (KrKh),/2T*/u*B*  (A  13) 

which  is  found  by  inspection  of  the 
normalized  equations.  This  characteristic  length  is  the 
scale  size  for  a  temperature  gradient  supported  by 
resistive  dissipation. 

There  is  an  additional  scale  size  for  variation 
due  to  the  change  in  magnetic  field  associated  with 
the  current  density  (that  drives  the  dissipation): 

dB  =  -|ij  (A  14) 

dx 

=  -  p  (E*  -  uB)  /  q 

In  terms  of  dimensionless  variables,  this  equation 
becomes: 

df  =  -  A  (1  +  1/Rm*  -  (Of )  03/2  (A  15) 
da 

where  an  additional  dimensionless  parameter 
is  obtained: 

A  =  n  (Kh  /  Kr  )l/2T*5/2  /  B*  (A16) 
=  Xc/(p*/nu*) 

that  relates  the  scale  size  for  thermal  conduction 
balancing  resistive  dissipation,  Xq,  to  that  for  which 
convection  balances  diffusion  of  magnetic  flux. 

Solution  of  the  set  of  normalized  equation  is 
obtained  by  integrating  in  the  upstream  direction 
from  the  starred  location,  where  conditions  are  taken 
as  f(0)  =  1,  9(0)  =  1,  and  T(0)  =  0,  (T  =  d9/da  ,  is  a 


dimensionless  temperature  gradient),  to  insure  that 
uniform  conditions  are  attained  in  the  limit  of  high 
magnetic  Reynolds  number.  The  actual  extent  of  the 
flow  field  is  not  prescribed,  but  is  determined  instead 
by  requirements  at  the  upstream  boundary  (e.g., 
necessary  heat  flux  to  establish  conditions  of  the 
entering  flow). 

For  the  one-dimensional  flow,  it  is  useful  to 
specify  a  sonic  condition  at  the  starred  location, 
rather  than  providing  values  for  the  mass  flow  that 
might  be  inconsistent  with  such  a  condition.  The 
necessary  value  of  u  is  then  given  (in  the  limit  of 
Rm  >:>l)by: 

u*2  =  yRT*  +  B*2  /  p*M  (A17) 

The  results  of  a  sample  calculation,  performed  using 
Mathematica,  are  displayed  in  Figures  A-l  to  A-3. 
Parameter  values  are  A  =  2.0,  and  P  =  0.1,  for  which 
P*  (=  2(y  -  1)PA  /  y  )  =  0.114.  The  local  value  of 
magnetic  Reynolds  number  is  (arbitrarily)  Rm*  =  10, 
and  the  specific  heat  ratio  is  y  -  1.4  .  Distances  are 
measured  upstream  of  the  sonic  point  by  the 
dimensionless  variable  b  =  x  /  Xc , 

To  return  to  dimensional  quantities,  it  is 
necessary  to  connect  the  results  of  the  normalized 
calculation  to  die  conditions  of  a  particular  thruster. 
For  example,  the  heat  flux  is: 

q  =  (KhT*7/2  /  Xc  )  05/2  T  (A  1 8) 

If  the  upstream  boundary  of  the  flow  corresponds  to 
the  entry  point  of  cold  propellant,  the  heat  flux  from 
the  discharge  must  be  sufficient  to  raise  the  total 
enthalpy  of  this  mass  to  the  initial  conditions  of  the 
discharge  flow.  For  purposes  of  illustration  here,  the 
necessary  heat  flux  may  be  written  as: 

qj  =  w  ( Q  +  CpT  i  +  uj^  /  2)  (A19) 

where  Q  is  the  chemical  energy  per  unit 
mass  (including  the  cost  of  vaporization,  dissociation 
and  ionization),  and  the  subscript  T  refers  to  the 
entry  station  of  the  flow.  The  characteristic 
temperature  T*  is  then  obtained  in  terms  of  the 
chemical  energy  per  unit  mass,  Q: 

t*  =  (  q  /  cp  y  [  ( e5^2  r  )i  /  p  -  01 

-  (y  - 1)(2  +  yP*)  to  i2  /  2yP*  ] 

(A20) 

In  Fig.  A-4,  the  denominator  of  Eqn.  21  is  displayed 
for  the  same  parameters  previously  used.  Note  that 
the  minimum  value  of  T*  corresponds  to  b  =  1.6  .  For 


a  Teflon  plasma,  fully  dissociated  into  singly-ionized 
constituents,  Q  is  about  62  eV/50  amu,  while  Cp 
would  be  21  eV/50  amu-eV  for  the  three  heavy- 
particles  and  three  electrons.  Thus,  Q/cp  —  3,  and  the 
minimum  value  of  T*  is  1.2  eV.  Higher  temperatures, 
however,  are  also  possible  and  would  be  chosen  in 
order  to  satisfy  other  conditions  of  the  thruster,  such 
as  the  operating  value  of  magnetic  field 

The  driving  source  for  the  thruster  can 
typically  be  characterized  in  terms  of  the  current 
supplied.  It  is  reasonable,  therefore,  to  attempt  to 
specify  thruster  operation  by  the  magnetic  field,  Bj, 
at  the  entrance  of  the  flow  field.  The  magnetic  field  at 
the  sonic  point  is  then  B*  =  Bj  /  f(bi)  .  The 
temperature  at  the  sonic  point  is  related  to  the 
magnetic  field  by: 

T*  -  {A  /  p(Kh  /  Kt)112}215  (Bj/fiCbi))275  (A21) 

Consistent  solution  requires  agreement  of  Eqns.  21 
and  22  for  a  specified  magnetic  field,  Bj.  In  the 
present  numerical  example,  this  occurs  at  bi  =  1.88. 
A  total  current  of  10  kA  over  a  2  cm  width  provides  a 
magnetic  field  of  Bi  =  0.63  Tesla  at  the  entrance 

implying  (with  f\  =  1.56)  a  value  of  B*  =  0.4  T  at  the 
sonic  point.  For  this  value,  the  characteristic 
temperature  may  be  found  in  terms  of  die  transport 
properties  of  the  plasma.  (For  an  ideal,  singly-ionized 
plasma,  the  values  of  Krand  Kh  are^: 

Kr  =  5.21  x  10'5  X  [W-m-eV3^2] 
and  Kh  =  7.46  x  104  /  A  [J  /  m-s-eV7/2] 

where  X  is  the  so-called  Coulomb  logarithm, 
and  temperatures  are  measured  in  eV.)  The 
characteristic  temperature  (with  X  =  10)  is  then: 

T*  =  {AB*/p  (Kh/Kr)1/2  (A22) 

=  7.8  eV 

With  the  magnetic  pressure,  and  plasma  temperature, 
the  mass  density  at  the  sonic  point  can  be  found  in 
terms  of  the  plasma-beta: 

p*  =  P*  (B*2  /  2p  )  /  RT*  (A23) 

=  8.3  x  10  kg/m3 

The  flow  speed  can  also  be  obtained  from 
the  sonic  condition  in  the  form: 

U*2  =  (B*2  /  p*p)  [  1  +  Yp*  /2  ]  (A24) 


For  the  numerical  values  previously  used,  u*  =  41 
km/sec.  The  characteristic  scale-size  is  0.93  mm,  so 
the  discharge  thickness  is  d  =  biXc  =  1.75  mm.The 

mass  flow  per  unit  area  is  p*u*  =  3.4  kg/m^-sec.  For 
a  cross-sectional  area  of  4  cm^,  the  mass  ejected  in 
10  ps  would  be  13.6  pg. 

Another  relationship  among  parameters: 

PA  =  jjp*CpT*  /  B*2  (A25) 

provides  the  speed  at  the  sonic  point  in  the  form: 


Fig.  A-l  Normalized  temperature,  0  vs  normalized 
distance,  b,  upstream  of  sonic  point 


Fig.  A-2  Normalized  magnetic  field,  f  vs  normalized 
distance,  b,  upstream  of  sonic  point 


u*=  Q1/2{[i+vp*/2]/[(e5/2r)I/p  -  0!  - 

(V  -  1)(2  +  Yp*)  CO  I2  /  2Yp*  ]PA  }  1/2 

(A26) 

which  displays  the  basic  scaling  with  Alfven 
critical  speed. 

While  the  numerical  results  in  the  present 
example  may  be  fortuitously  close  to  values  observed 
in  various  PPT  experiments,  accurate  predictions 
require  modeling  based  on  the  actual  behavior  of  the 
propellant  in  the  hill  two-  (and  three-)  dimensional, 
unsteady  environment  of  the  thruster. 


Fig.  A-3  Normalized  flow  speed,  a)  vs  normalized 
distance,  b,  upstream  of  sonic  point 


Fig.  A-4:  [(  05/2  T  )i  /  P  -  0 1  -  (Y  -  1)(2  +  Yp*)  w  i2 
/  2Yp*]  of  Eqn.  20  vs  normalized 
distance,  b,  upstream  of  sonic  point 
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ABSTRACT 

The  chemical  composition  of  tetrafluoroethylene 
(C2F4)  is  calculated  with  a  two-temperature  LTE 
formulation.  Twenty-five  chemical  species  are 
included  in  the  analysis.  The  equilibrium  constants 
are  calculated  using  the  most  recent  spectroscopic 
data  available.  Calculations  are  performed  for 
pressures  from  0.001  atm  to  1.0  atm  and  for 
temperature  ranges  of  0.05  ev  to  10  ev  for  both 
heavy  particle  and  electron  temperatures. 


INTRODUCTION 

Knowledge  of  the  chemical,  thermodynamic,  and 
transport  properties  of  a  gas  is  required  in  almost  any 
gasdynamic  analysis.  Accurate  thermochemical  and 
transport  properties  become  particularly  important  in 
high-temperature  applications  such  as  the  pulsed 
plasma  thruster.  In  this  paper,  we  shall  concentrate 
on  calculating  the  equilibrium  composition  of  a  gas 
mixture.  This  is  the  necessary  first  step  for 
determining  the  thermodynamic  and  transport 
properties  of  a  gas. 

There  are  three  primary  thermochemical  states 
possible  for  a  gas.  A  calorically  perfect  gas  has 
specific  heats  that  are  constant,  and  the  enthalpy  and 
internal  energy  are  only  functions  of  temperature.  A 
thermally-perfect  gas,  in  which  variable  vibrational 
and  electronic  excitation  are  taken  into  account,  has 
specific  heats,  enthalpy,  and  internal  energy  that  are 
all  functions  of  temperature.  If  the  conditions  are 
right  for  chemical  reactions  to  occur,  then  we  can 
treat  the  gas  as  an  equilibrium  chemically-reacting  gas 


for  which  properties  area  all  functions  of  temperature 
and  pressure.  Even  this  can  be  generalized  by  stating 
that  the  reacting  gas  is  in  local  thermodynamic 
equilibrium  (LTE).  This  means  that  a  local 
Boltzmann  distribution  exists  at  each  point  in  the 
flow  at  the  local  temperature.  We  will  extend  this 
statement  further  for  the  case  of  a  two  temperature 
LTE  gas  modeled  here.  In  this  paper  we  calculate  the 
chemical  composition  of  the  gaseous  Teflon  monomer 

(C2F4). 


POSSIBLE  SPECIES,  REACTIONS,  AND 
EQUILIBRIUM  EQUATIONS 

In  this  paper,  we  calculate  the  chemical 
composition  of  tetrafluoroethylene  (C2F4).  The 
analysis  will  include  vibrational  and  electronic 
excitation,  dissociation,  first  molecular  ionization, 
and  first  through  fourth  monatomic  ionization. 
Throughout  of  the  analysis,  we  shall  assume  a  perfect 
gas,  where  intermolecular  forces  are  non-existent  or 
negligible.  This  might  seem  a  strange  assumption 
when  the  gas  is  in  the  plasma  state  due  to  the 
presence  of  Coulomb  collisions,  but  it  is  a  widely 
used  and  accepted  approximation.1 

For  a  polyatomic  base  gas,  C2F4  in  our  case,  with 
the  possibility  of  undergoing  full  dissociation, 
singular  molecular  ionization,  and  up  to  fourth 
monatomic  ionization,  we  first  assume  there  are 
twenty-five  possible  chemical  species,  which  are 
C2F4,  c2f2,  cf2,  cf2\  cf3,  CF3\  CF4,  C2,  CF, 
CF+,  F2,  F2\  Cz  (Z— 1,4),  Fz  (Z=-l,4),  and  e  's. 
For  a  gas  containing  twenty-five  chemical  species, 
which  is  composed  of  three  elements  (C,  F,  e  ),  we 
are  required  to  have  twenty-two  (25-3=22) 
independent  chemical  reaction  equations  (laws  of 
mass  action).  The  reactions  considered  here  are 


CF2  ~CF  +  F 

(1) 

CF  ~C  +  F 

(2) 

F2  ~2F 

(3) 
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CF2  ~  CFl  +  e~ 

(4) 

CF  ~  CF'  +  e~ 

(5) 

4 

z~l  Cz  +  e~  1 

(6) 

2-0 

C2F4  -  2  CF2 

(7) 

f2  *  f;  *  e- 

(8) 

formulations,  they  may  also  be  put  in  terms  of 
concentrations,  K,.,  or  number  densities,  Kn.  It  is 
important  to  note  in  the  above  equations  that  the 
equilibrium  constants  are  written  as  functions  of 
temperature  only,  as  most  authors  point  out. 
However,  they  may  be  functions  of  two  or  more  state 
variables  depending  on  whether  such  things  as 
thermal  non-equilibrium  assumptions  or  electronic 
partition  function  cutoff  is  taken  into  account.2 

In  addition  to  twelve  independent  equations 
relating  the  twenty-five  unknown  partial  pressures, 
we  need  three  more  equations  to  solve  for  the  gas 
composition.  The  three  chosen  are;  conservation  of 
nuclei,  Dalton's  Law,  and  charge  neutrality.  The 
ideal  thermal  gas  law  for  each  species  is  written  in 
the  form 


Pt  -  nikTi  <15> 


Fzl~Fz  +  e-  | 

2=0 

(9) 

where 

P 

-E  Pi 

(16) 

C2Fa  *■*  c2f2  +  f2 

(10) 

For  charge  neutrality,  we  have 

C2Fa  ~  CF3  +  CF 

(ID 

25 

£*»< 
1  =  1 

2  =  0 

(17) 

C2Fa  **  CF4  +  c 

(12) 

In  terms  of  partial  pressures,  this  becomes  (  for  ideal 
gases  ) 

C2  •*>  2C 

03) 

25 

E  z 

i=l 

1] 

T, 

l) 

Pi2  =  0 

08) 

In  actuality,  there  are  other  possible  reactions  that 
could  yield  the  same  chemical  species.  But,  for  an 
equilibrium  calculation,  the  reactions  chosen  are 
arbitrary  as  long  as  they  are  linearly  independent  and 
account  for  all  possible  species. 

Writing  these  reactions  in  terms  of  equilibrium 
relations  for  the  partial  pressures,  we  have 

kjt)  -  n  p'‘  <i4> 


For  conservation  of  nuclei,  we  write 


where 


25 


("c)0  -  E 

i- 1 


CN 


(20) 


where  the  Kpj  are  the  equilibrium  constants  for  the 
reaction  (j)  at  the  equilibrium  temperature  T,  in  terms 
of  the  partial  pressures.  Using  the  appropriate 


and 


(H  -  E  Nx„  <21> 

where  (nC2Fi)o  is  the  total  number  of 
tetrafluoroethylene  molecules  available  for 
dissociation  and  ionization  (ie.  the  number  of  CjF* 
molecules  present  if  the  gas  was  non-reacting  at  some 
initially  low  temperature).  Dividing  Eq.  (20)  by  Eq. 
(21)  and  utilizing  Eq.  (19)  gives  us  the  nuclei 
conservation  statement,  where  the  number  densities 
are  related  to  the  partial  pressures  by  Eq.  (15). 


CALCULATION  OF  THE  EQUILIBRIUM 
CONSTANT  -  PARTITION  FUNCTIONS 

To  solve  the  system  of  equations,  we  only  need 
values  for  the  equilibrium  constants  which  may  be 
calculated  from  equilibrium  statistical  mechanics.  In 
terms  of  partition  functions  Qit  the  law  of  mass  action 
for  a  general  system  is 

-Ae,  (22) 

k/t,  =  IK  - «  "  n<?." 


or  alternatively,  substituting  ^  =  Nj/V  we  have 


(23) 


where  vs  is  the  stoichiometric  mole  number  for 
species  (i),  that  is,  the  coefficients  in  the  balanced 
chemical  equation,  Ae0  is  the  reaction  energy 
(change  in  zero-point  energy)  and  Qj  is  the  total 
partition  function  for  species  (i).  Thus,  for  a  given 
reaction  and  thermodynamic  state,  the  only  unknowns 
in  Eq.  (23)  are  the  Q/s. 

For  a  system  in  thermodynamic  equilibrium,  we 
have 


which  gives  the  number  of  particles  Nj*  in  energy 
level  Cj  with  gj  degenerate  states.  We  define  the 
partition  function,  Q,  as  the  sum  in  the  denominator 
of  Eq.  (24)  which  is,  in  general,  a  function  of  T  and 


Q  -=  E  s/tr 


(25) 


V.  It  is  typical  to  express  the  total  energy  as  the  stun 
of  translational  and  internal  energies.  Note  that  Eqs. 
(23)  and  (24)  contain  only  one  temperature.  For  the 
two  temperature  case  considered  in  this  research  we 
make  the  assumption  that  the  heavy-particle  gas 
composed  of  neutrals  and  ions,  has  a  Maxwellian 
distribution  in  velocities  and  a  Boltzman  distribution 
in  energies  at  a  heavy-particle  temperature  ,  T.  The 
electron  gas,  composed  of  both  free  and  bound 
electrons  is  in  equilibrium  with  an  electron 
temperature  Te  defined  by  their  Maxwellian  velocity 
distribution.  Note  that,  in  this  analysis,  we  are 
ignoring  the  interaction  between  electronic  and 
vibrational  states.  Thus,  we  have  defined  a  two- 
temperature  LTE  situation.  For  a  molecule  we  have 

e=c  +  e  .  +  e  .,  +  e  ,  (26) 

*  trans  rot  vtb  el 


and  for  an  atom 


e 


e  +  e  , 

trans  el 


(27) 


where  c  is  the  sensible  energy,  measured  above  the 
zero-point  energy  e0.  Quantum  mechanics  gives  us 
theoretical  values  for  the  quantized  energies  of  a 
particle,  at  least  for  the  translational,  rotational,  and 
vibrational  modes.3  Along  with  the  assertion  that 
particle  energy  is  simply  the  sum  of  the  modal 
energies,  that  is,  the  internal  energies  are  uncoupled, 
which  is  a  consequence  of  the  more  fundamental 
assumption  of  a  separable  Hamiltonian,  the  partition 
function  is  expressed  as  the  product  of  the  modal 
partitions  Qj,  where 

q  -  n  <?,  <28> 


(24) 


with  j  extending  over  all  modes.  Armed  with  the 
quantized  values  for  the  modal  energies,  and  the 
associated  degeneracies  we  can  calculate  the  modal 
partition  functions  which  are  given  here  in  reduced 
form,  without  proof  as4 


Qtrans 


InmkT 


h2 


xV 


(29) 


of  ionization  potential  lowering  is  the  subject  of  some 
controversy  and  was  explored  in  detail  in  a  previous 
work.2  Results  from  that  work  give  the  correct  cutoff 
criterion  as 


■rot 


ZjflkT 

oh2 


(30) 
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cutoff  _ 


1  - 


,(33) 


X  n„ 


Qvib  ~ 
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(31) 


where  a  is  a  factor  which  arises  from  the  symmetry 
requirements  of  the  wave  function  in  the  exchange  of 
an  identical  particle.  It  is  equal  to  1  for  heteronuclear 
molecules  (ex.  CF),  and  equal  to  2  for  homonuclear 
molecules  (ex.  F2).  For  electronic  energy  there  is  no 
closed  form  general  expression  for  the  quantized 
energy  levels.  Thus  the  electronic  partition  function 
must  be  left  as  an  infinite  series  in  the  form 

=  (32) 
j-0 


where  a=Zcff*e2/2*IP.  The  lowered  ionization 
potential  is  given  by 


IP  =  IP 


1  - 


n 


cutoff 


(34) 


The  Kn’s  are  converted  to  Kp’s  using  the  relation 

kjt) = <*Nn T:' 

The  required  molecular  and  atomic  data,  which  are 
too  numerous  to  give  here,  are  taken  from  the  works 
of  Chase,7  Moore, 8,9  Rosenstock,10  Herzberg,11 
Buckely,12  and  Paulino  and  Squires.13 


,(35) 


Equation  (32)  is  the  true  theoretical  representation  of 
the  electronic  partition  function  for  an  isolated 
particle.  In  theory  there  are  in  infinite  number  of 
electronic  levels  extending  from  the  ground  state 
energy  (e0  =  0)  to  the  ionization  potential,  which  is 
the  amount  of  energy  needed  to  remove  an  electron 
from  its  ground  state  to  infinity  (bound-free 
transition).  The  electronic  partition  function  is  a 
diverging  series  because  although  the  energy 
approaches  a  finite  limit,  the  degeneracy  increases  as 
the  square  of  the  principal  quantum  number,  so  the 
series  diverges.5  For  any  general  polyatomic  molecule 
of  N  atoms,  if  we  still  assume  a  separable 
Hamiltonian  then  we  can  factor  the  partition  function 
as  in  Eq.  (28)  with  the  product  extending  not  only 
over  all  fundamental  modes  but  also  over  all  modal 
degrees  of  freedom.6 

In  actuality,  the  electronic  series  in  not  infinite 
because  a  particle  in  the  real  world  is  never  truly 
isolated.  Due  to  various  interparticle  interactions  that 
arise  in  any  finite  density  medium,  the  series  will 
actually  terminate  at  some  principal  quantum  number, 
jjcwoff  -pjje  evaluation  of  ncumff  and  its  associated  effect 


SOLUTION  OF  THE  EQUILIBRIUM 
EQUATIONS 

The  equations  given  in  the  previous  section  to 
calculate  the  equilibrium  composition  in  terms  of 
partial  pressures,  given  T  ,  Te  and  P,  provide  a 
closed  set  of  twenty-five  coupled  nonlinear  algebraic 
equations  for  which  there  is  no  analytic  solution  and 
numerical  methods  must  be  used.  The  numerical 
solution  of  systems  of  nonlinear  equations  is 
universally  very  difficult  and  is  a  topic  of  current 
research.  As  of  yet,  there  are  no  appropriate 
numerical  methods  for  solving  coupled  nonlinear 
systems  of  algebraic  equations  from  arbitrary  starting 
vectors.14  One  of  the  most  commonly  used  methods 
and  the  one  used  previously  for  nitrogen  research,2  is 
the  Newton-Raphson  method.  The  Newton-Raphson 
method  usually  exhibits  excellent  convergence 
qualities  when  the  starting  vector  (initial  guesses  for 
die  roots)  is  near  the  actual  root.  For  a  homonuclear 
diatomic  gas  a  good  initial  guess  can  usually  be 
obtained  through  Saha-type  statistical  arguments  using 


weighted  averages  due  to  limited  reactive 
simultaneity.  In  this  case  though,  since  we  are 
dealing  with  a  polyatomic  heteronuclear  gas,  it  was 
initially  assumed  that  we  could  not  get  close  enough 
to  the  root  for  Newton-Raphson  to  work  properly. 
Many  hybrid  techniques  have  been  proposed  but  most 
fail  when  the  Jacobian  becomes  singular,  or  at 
stationary  point.1516  The  solution  procedure  initially 
chosen  for  this  research  is  the  one  proposed  by 
Powell.17  Powell’s  technique  exhibits  almost 
guaranteed  convergence  even  for  poor  initial  guesses. 
It  also  has  the  ability  to  correctly  handle  stationary 
points.  The  price  of  this  behavior  is  that  the 
convergence  is  linear  until  very  close  to  the  root  then 
it  converges  quadradcally  like  the  standard  Newton- 
Raphson.  Thus,  the  total  number  of  iterations 
required  is  greatly  increased  resulting  in  much  longer 
runtimes,  up  to  200  minutes  per  isobar  on  a  Sun 
Supersparc  .  Thus,  it  was  decided  that  the  Newton- 
Raphson  method  would  be  used  combined  with  a 
raster  processing  iteration  procedure. 

Given  T,  Te,  and  P,  an  initial  equilibrium 
composition  is  approximated  using  simple  Saha-type 
arguments  assuming  totally  uncoupled  physical 
process,  then  partial  pressures  are  backed  out  using 
weighted  approximations  or  the  initial  guess  was  set 
equal  to  die  previous  converged  root  at  the  last 
temperature  and  pressure  and  stepped  up  in 
temperature  on  each  isobar  in  small  increments  of 
from  0. 1  to  10  K  for  heavy  particle  temperatures  and 
1  K  to  100  K  for  electron  temperatures.  The 
completely  closed  system  is  solved  at  each  point 
using  the  Newton-Raphson  procedure,  which  iterates 
until  the  sum  of  the  absolute  values  of  the  corrections 
is  less  than  the  chosen  tolerance  of 


eP 


(36) 


where  e  =  10-4  or  until  20,000  iterations  are 
performed.  Convergence  to  the  tolerance  is  usually 
obtained  quickly.  An  additional  convergence  check  is 
performed  after  each  completed  iteration  so  that  the 
total  isobaric  deviation  is  not  allowed  to  exceed  5%. 
The  correct  composition  is  now  known  to  within 
desired  accuracy.  The  mixture  molecular  is  was 
computed  using  the  expression 

25 

MMG  =  Y,  XiMMi  (37) 

i=  1 


where  X5  is  the  mole  fraction  of  species  i,  obtained 
using 


(38) 


where  n^,  is,  of  course,  just  the  total  number  of 
particles. 


RESULTS  AND  DISCUSSION 

The  chemical  composition  in  terms  of  the  mole 
fractions  at  1  atm.  and  Te/T  =  1.0  from  300  to 
30,000  K  is  shown  in  Fig.  1.  We  see  the  same  rapid 
dissociation  of  CyF4  as  observed  by  Paulino  and 
Squires13  due  to  the  inherent  weakness  of  the  carbon 
double  bond  in  this  molecule.  The  results  are 
consistent  with  those  obtained  by  Kovitya.18  Figure  2 
shows  a  better  view  of  the  composition  in  the  low 
temperature  region  where  molecules  dominate.  In  our 
case  we  see  that  at  about  800  K  CjF4  partially 
dissociates  into  QF2  and  CF4,  almost  completely 
recombines  as  temperature  increases,  then  dissociates 
into  CF2.  Dissociation  of  CF2  begins  to  at  about 
4,000  K  and  is  almost  completely  dissociated  by 
7,500  K.  The  primary  dissociation  products  are  C 
and  F  which  reach  their  maxima  at  about  7,500  K. 
Past  this  temperature  ionization  begins  to  occur  and 
singly  ionized  C  and  F  and  electrons  dominate  the 
composition.  Second  ionization  begins  to  occur 
between  22,500  and  25,000  K.  Figures  3  through  8 
show  representative  species  compositions  of  C,  C+, 
F,  F+,  e‘  and  CF2  at  two  different  pressures,  1  atm 
and  .  1  atm,  for  the  case  of  thermal 
nonequilibrium.They  are  given  for  four  heavy 
isotherms  of  .05,  .1,  1,  and  10  ev.  Here  we  see 
drastically  different  behavior  due  to  the  highly 
energetic  electrons  at  their  elevated  temperature. 
Dissociation  of  CF2  happens  rapidly,  other  molecules 
exhibited  similar  behavior.  Single  and  multiple 
ionization  of  C  and  F  follows  the  expected  pattern  of 
following  the  electron  temperature,  consistent  with 
our  original  formulation.  Figure  7,  the  electron 
partial  pressure,  shows  we  reach  regions  where 
dPe/dTe  goes  to  zero.  The  composition  is  constant 
after  this  for  we  had  assumed  that  the  maximum  ionic 
charge  state  was  four.  This  suggests  that  for  a  correct 
chemical  model  at  elevated  electron  temperatures  we 
may  need  to  extend  the  possible  ionic  species  to  a 
greater  charge  value.  The  molecular  ionic  and 
electronegative  species  were  found  to  exist  in  only 


very  small  amounts  at  all  Te/T  values.  In  future  work 
it  is  planned  to  include  multi-phase  species  in  like 
amorphous  carbon.  Molecular  effects  that  are  being 
evaluated  for  inclusion  in  future  models  include 
anharmonic  vibrations  and  internal  energy  mode 
coupling.  Also  research  is  being  performed  to 
develop  the  complete  reactive  thermodynamic  and 
transport  models. 
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Figure  2.  Composition  of  Teflon  for  Te/T  =  1.0  at 
1  atm.  Lower  temperature  range. 


Figure  1 .  Composition  of  Teflon  for 
Te/T  =  1.0  at  1  atm. 


Figure  3.  Partial  Pressure  of  C  vs.  Te. 
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Partial  Pressure  of  C+  vs.  Te 


Partial  Pressure  of  F+  vs.  Te 


Figure  4.  Partial  Pressure  of  C+  vs.  Te. 


Figure  6.  Partial  Pressure  of  F+  vs.  Te. 


Figure  5.  Partial  Pressure  of  F  vs.  Te. 


Figure  7.  Partial  Pressure  of  electrons  vs.  Te. 
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Figure  8.  Partial  Pressure  of  CF2  vs.  Te. 
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ABSTRACT 

The  equilibrium  chemical  composition,  thermodynamic,  and  transport  properties  of  nitrogen  is 
calculated  over  a  wide  range  of  thermodynamic  state  conditions.  Six  chemical  species  are 
included  in  the  analysis.  The  equilibrium  constants  are  calculated  using  the  most  recent 
spectroscopic  data  available.  Calculations  are  performed  for  pressures  from  0.001  atm  to  1 ,000.0 
atm  and  for  temperature  ranges  of  300  K  to  30,000  K. 


INTRODUCTION 

Knowledge  of  the  chemical,  thermodynamic,  and  transport  properties  of  a  gas  is  required  in 
almost  any  gasdynamic  analysis.  Accurate  thermochemical  and  transport  properties  become 
particularly  important  in  high-temperature  applications  such  as  the  pulsed  plasma  thruster  and 
other  high  temperature  devices.  In  this  paper,  we  shall  concentrate  on  calculating  the  equilibrium 
composition,  thermodynamic,  and  transport  properties  of  a  gas  mixture,  in  this  case  nitrogen. 

There  are  three  primary  thermochemical  states  possible  for  a  gas.  A  calorically  perfect  gas  has 
specific  heats  that  are  constant,  and  the  enthalpy  and  internal  energy  are  only  functions  of 
temperature.  A  thermally-perfect  gas,  in  which  variable  vibrational  and  electronic  excitation  are 
taken  into  account,  has  specific  heats,  enthalpy,  and  internal  energy  that  are  all  functions  of 
temperature.  If  the  conditions  are  right  for  chemical  reactions  to  occur,  then  we  can  treat  the  gas 
as  an  equilibrium  chemically-reacting  gas  for  which  properties  area  all  functions  of  temperature 
and  pressure.  Even  this  can  be  generalized  by  stating  that  the  reacting  gas  is  in  local 
thermodynamic  equilibrium  (LTE).  This  means  that  a  local  Boltzmann  distribution  exists  at  each 
point  in  the  flow  at  the  local  temperature.  Recent  research,  both  in  the  areas  of  advanced 
propulsion  and  high-temperature  gas  dynamics  has  demonstrated  the  need  for  accurate  data  on 
the  chemical,  thermodynamic,  and  transport  properties  of  such  reacting  gases.  In  this  paper  we 
analyze  the  thermophysical  properties  of  nitrogen  with  the  inclusion  of  translational,  rotational, 
vibrational,  and  electronic  excitation. 


POSSIBLE  SPECIES,  REACTIONS,  AND 
EQUILIBRIUM  EQUATIONS 

In  this  section,  we  calculate  the  chemical  composition  of  nitrogen.  The  analysis  includes 
vibrational  and  electronic  excitation,  dissociation,  first  molecular  ionization,  and  first  through 
second  monatomic  ionization.  Throughout  of  the  analysis,  we  assume  a  perfect  gas,  where 


intermolecular  forces  are  non-existent  or  negligible.  This  might  seem  a  strange  assumption  when 
the  gas  is  in  the  plasma  state  due  to  the  presence  of  Coulomb  collisions,  but  it  is  a  widely  used 
and  accepted  approximation.1 

For  a  diatomic  base  gas,  N2  in  our  case,  with  the  possibility  of  undergoing  full  dissociation, 
singular  molecular  ionization,  and  up  to  second  monatomic  ionization,  we  first  assume  there  are 
six  possible  chemical  species,  which  are  N2,  N2+,  Nz  (z=0,2)  ,  and  e'  's.  For  a  gas  containing 
six  chemical  species,  which  is  composed  of  two  elements  (N,  e‘),  we  are  required  to  have  four 
(6-2=4)  independent  chemical  reaction  equations  (laws  of  mass  action).  The  reactions  considered 
here  are 


N2  +  e"2  ~2N 
N2  +  ef2  N;  +  e- 
N  +  ef^  N *  +  e~ 


(1) 

(2) 

(3) 


AT  +  N++  +  e~ 

In  actuality,  there  are  other  possible  reactions  that  could  yield  the  same  chemical  species.  But, 
for  an  equilibrium  calculation,  the  reactions  chosen  are  arbitrary  as  long  as  they  are  linearly 
independent  and  account  for  all  possible  species. 

Writing  these  reactions  in  terms  of  equilibrium  relations  for  the  partial  pressures,  we  have 

kjt>  -  n  pi'  <s> 


where  the  KpJ  are  the  equilibrium  constants  for  the  reaction  (j)  at  the  equilibrium  temperature 
T,  in  terms  of  the  partial  pressures.  Using  the  appropriate  formulations,  they  may  also  be  put 
in  terms  of  concentrations,  Kc,  or  number  densities,  Kn. 

In  addition  to  four  independent  equations  relating  the  six  unknown  partial  pressures,  we  need 
two  more  equations  to  solve  for  the  gas  composition.  The  two  chosen  are  Dalton's  Law  and 
charge  neutrality.  The  ideal  thermal  gas  law  for  each  species  is  written  in  the  form 

p.  =  nJkT  (6) 

where 

p  =  £  Pl  (7) 


For  charge  neutrality,  we  have 
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£zv  =  o 

i=l 

In  terms  of  partial  pressures,  this  becomes  (  for  ideal  gases  ) 

E  zp,‘  - 0 

i=l 


(9) 


CALCULATION  OF  THE  EQUILIBRIUM 
CONSTANT  -  PARTITION  FUNCTIONS 


To  solve  the  system  of  equations,  we  only  need  values  for  the  equilibrium  constants  which 
may  be  calculated  from  equilibrium  statistical  mechanics.  In  terms  of  partition  functions  Q;,  the 
law  of  mass  action  for  a  general  system  is 


"A  £a 
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or  alternatively,  substituting  n,  =  N/V  we  have 


Kn(D 


(ID 


where  Vj  is  the  stoichiometric  mole  number  for  species  (i),  that  is,  the  coefficients  in  the 
balanced  chemical  equation,  Ae0  is  the  reaction  energy  (change  in  zero-point  energy)  and  Q,  is 
the  total  partition  function  for  species  (i).  Thus,  for  a  given  reaction  and  thermodynamic  state, 
the  only  unknowns  in  Eq.  (11)  are  the  Q,'s. 

For  a  system  in  thermodynamic  equilibrium,  we  have 

NJ  =  N  R'C  -  02) 

IL 
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which  gives  the  number  of  particles  N/  in  energy  level  Cj  with  gj 
the  partition  function,  Q,  as  the  sum  in  the  denominator  of  Eq. 

Q  =  £  Sje" 


degenerate  states.  We  define 
(12)  which  is,  in  general,  a 

(13) 


function  of  T  and  V.  It  is  typical  to  express  the  total  energy  as  the  sum  of  translational  and 


internal  energies.  Note  that,  in  this  analysis,  we  are  ignoring  the  interaction  between  electronic 
and  vibrational  states.  Thus,  we  have  defined  a  one-temperature  LTE  situation.  For  a  molecule 
we  have 
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and  for  an  atom 
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where  e  is  the  sensible  energy,  measured  above  the  zero-point  energy  e0.  Quantum  mechanics 
gives  us  theoretical  values  for  the  quantized  energies  of  a  particle,  at  least  for  the  translational, 
rotational,  and  vibrational  modes.2  Along  with  the  assumption  that  particle  energy  is  simply  the 
sum  of  the  modal  energies,  that  is,  the  internal  energies  are  uncoupled,  which  is  a  consequence 
of  the  more  fundamental  assumption  of  a  separable  Hamiltonian,  the  partition  function  is 
expressed  as  the  product  of  the  modal  partitions  Qj,  where 

<?  =  !!<?;  (16) 


with  j  extending  over  all  modes.  Armed  with  the  quantized  values  for  the  modal  energies,  and 
the  associated  degeneracies  we  can  calculate  the  modal  partition  functions  which  are  given  here 
in  reduced  form,  without  proof  as3 
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where  a  is  a  factor  which  arises  from  the  symmetry  requirements  of  the  wave  function  in  the 
exchange  of  an  identical  particle.  It  is  equal  to  1  for  heteronuclear  molecules  (ex.  NO),  and  equal 
to  2  for  homonuclear  molecules  (ex.  N2).  For  electronic  energy  there  is  no  closed  form  general 
expression  for  the  quantized  energy  levels.  Thus  the  electronic  partition  function  must  be  left  as 
an  infinite  series  in  the  form 
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Eq.  (20)  is  the  true  theoretical  representation  of  the  electronic  partition  function  for  an  isolated 
particle.  In  theory  there  are  in  infinite  number  of  electronic  levels  extending  from  the  ground 
state  energy  (e0  =  0)  to  the  ionization  potential,  which  is  the  amount  of  energy  needed  to 
remove  an  electron  from  its  ground  state  to  infinity  (bound-free  transition).  The  electronic 
partition  function  is  a  diverging  series  because  although  the  energy  approaches  a  finite  limit,  the 
degeneracy  increases  as  the  square  of  the  principal  quantum  number,  so  the  series  diverges.4,5 

In  actuality,  the  electronic  series  in  not  infinite  because  a  particle  in  the  real  world  is  never 
truly  isolated.  Due  to  various  interparticle  interactions  that  arise  in  any  finite  density  medium, 
the  series  will  actually  terminate  at  some  principal  quantum  number,  n0^.  The  evaluation  of 
ncmoff  jjjyj  jts  associated  effect  of  ionization  potential  lowering  is  the  subject  of  some  controversy 
and  was  explored  in  detail  in  a  previous  work.6  Results  from  that  work  give  the  correct  cutoff 
criterion  as 
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where  a=Zeff*e2/2*IP.  The  lowered  ionization  potential  is  given  by 
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The  Kn's  are  converted  to  K's  using  the  relation 
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(23) 


The  required  molecular  and  atomic  data,  which  are  too  numerous  to  give  here,  are  taken  from 
the  works  of  Herzberg7  and  Moore8. 


SOLUTION  OF  THE  EQUILIBRIUM 
EQUATIONS 

The  equations  given  in  the  previous  section  to  calculate  the  equilibrium  composition  in  terms 
of  partial  pressures,  given  T  ,  Te  and  P,  provide  a  closed  set  of  six  coupled  nonlinear  algebraic 
equations  for  which  there  is  no  analytic  solution  and  numerical  methods  must  be  used.  The 
numerical  solution  of  systems  of  nonlinear  equations  is  universally  very  difficult  and  is  a  topic 
of  current  research.  As  of  yet,  there  are  no  appropriate  numerical  methods  for  solving  coupled 
nonlinear  systems  of  algebraic  equations  from  arbitrary  starting  vectors.14  One  of  the  most 


commonly  used  methods  is  the  Newton-Raphson  method.  The  Newton-Raphson  method  usually 
exhibits  excellent  convergence  qualities  when  the  starting  vector  (initial  guesses  for  the  roots)  is 
near  the  actual  root.  For  a  homonuclear  diatomic  gas  a  good  initial  guess  can  usually  be  obtained 
through  Saha-type  statistical  arguments  using  weighted  averages  due  to  limited  reactive 
simultaneity.  Another  technique,  known  as  Powell's  technique  exhibits  almost  guaranteed 
convergence  even  for  poor  initial  guesses.10  It  also  has  the  ability  to  correctly  handle  stationary 
points.  The  price  of  this  behavior  is  that  the  convergence  is  linear  until  very  close  to  the  root 
then  it  converges  quadratically  like  the  standard  Newton-Raphson.  Thus,  the  total  number  of 
iterations  required  is  greatly  increased  resulting  in  much  longer  runtimes.  Thus,  it  was  decided 
that  the  Newton-Raphson  method  would  be  used  combined  with  a  raster  processing  iteration 
procedure. 

Given  T  and  P,  an  initial  equilibrium  composition  is  approximated  using  simple  Saha-type 
arguments  assuming  totally  uncoupled  physical  process,  then  partial  pressures  are  backed  out 
using  weighted  approximations  or  the  initial  guess  was  set  equal  to  the  previous  converged  root 
at  the  last  temperature  and  pressure  and  stepped  up  in  temperature  on  each  isobar  in  small 
increments  of  from  100  to  300  K  for  heavy  particle  temperatures.  The  completely  closed  system 
is  solved  at  each  point  using  the  Newton-Raphson  procedure,  which  iterates  until  the  sum  of  the 
absolute  values  of  the  corrections  is  less  than  the  chosen  tolerance  of 


where  e  =  10"4  or  until  20,000  iterations  are  performed.  Convergence  to  the  tolerance  is  usually 
obtained  quickly.  An  additional  convergence  check  is  performed  after  each  completed  iteration 
so  that  the  total  isobaric  deviation  is  not  allowed  to  exceed  5%.  The  correct  composition  is  now 
known  to  within  desired  accuracy.  The  mixture  molecular  is  was  computed  using  the  expression 
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MMG  =  Y,X<  MMi  (25) 

i  =  1 


where  X;  is  the  mole  fraction  of  species  i,  obtained  using 


(26) 


where  ntot  is,  of  course,  just  the  total  number  of  particles. 


CALCULATION  OF  THE  THERMODYNAMIC  PROPERTIES 


For  a  given  temperature  and  pressure,  the  complete  gas  composition  can  be  calculated  using 
the  techniques  in  the  previous  section.  Often,  however,  it  is  necessary  to  have  the  thermodynamic 
propertied  of  the  gas.  In  this  research,  the  calculated  thermodynamic  properties  are  specific 
enthalpy  (h),  specific  internal  energy  (e),  the  specific  heats  (cp  and  Cv),  and  consequently  the  ratio 


of  specific  heats  (y). 


a.  Calculation  of  the  specific  internal  energy  (e)  and  specific  enthalpy  (h) 

On  a  microscopic  scale,  a  particle  (atom  or  molecule)  may  possess  energy  due  to  its 
translational,  rotational,  vibrational,  and  electronic  modes.  This  is  expressed  as 
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where  e  is  the  total  specific  internal  energy.  Note  that  e  is  actually  the  total  specific  sensible 
internal  energy  and  is  measured  above  the  zero  point.  In  terms  of  partition  functions  the  energy 
and  enthalpy  is  expressed  as3 

e  =  (28) 


ft  =  RT  +  <29> 


with  the  relation 

h  =  e  +  Pv  (30) 


holding  for  any  type  of  gas,  whether  calorically  perfect,  thermally  perfect,  or  chemically 
reacting.  The  expressions  for  the  modal  partition  functions  are  given  in  the  previous  section. 
Taking  the  derivatives  of  the  logarithms  and  using  Eq.(28)  we  get 
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Again,  for  the  electronic  energy,  the  lack  of  a  closed  form  partition  function  leaves  us  with 
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where  R;  =  RUGC/MMi  is  the  specific  gas  constant.  The  translational  and  rotational  modes  are 
assumed  to  be  fully  excited. 

For  a  mixture  of  N  species,  the  total  mixture  sensible  internal  energy  is 
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where  c,  is  the  mass  fraction  of  species  i,  given  by 

MMi  P .  MM. 
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Similarly,  the  specific  enthalpy  for  a  mixture  of  N  species  is  given  by 


N 


h  =  £  cihi 
»=i 


(37) 


K  =  +  RJ  +  (A  hf)\ 


(38) 


where  again,  c,  is  the  mass  fraction  and  hj  is  the  specific  enthalpy  for  species  i.  The  internal 
energy  at  each  required  temperature  and  pressure  is  calculated  using  Eqs.  (31)-(34)  combined 
with  Eqs.  (35)  and  (36). The  solution  procedure  required  a  previous  call  to  the  equilibrium 
composition  subroutine  to  input  the  partial  pressures  and  mixture  molecular  mass.  The  electronic 
derivatives  appearing  in  Eqs.  (31)  through  (34)  are  calculated  using  a  4th-order  four-point 
Richardson  extrapolation  with  AT  =  .OIK.  The  derivatives  are  taken  with  a  frozen  principal 
quantum  number  cutoff.  This  was  done  for  ease  of  calculation  and  to  avoid  the  problem  other 
authors4,5  experienced  with  large  computational  derivatives  being  obtained  if  the  temperature  and 
pressure  state  corresponded  to  a  partition  function  jump  due  to  the  addition  or  subtraction  of 
another  electronic  level.  The  specific  enthalpy  is  not  computed  using  Eq.(37)  ,  rather  Eq.(30) 
is  used  written  in  the  form 


(39) 


The  mass  fractions  are  input  to  the  thermodynamic  subroutine,  thus  are  known.  The  (dh;/dT)p 
for  each  species  is  calculated  using  Eq.(38)  together  with  the  appropriate  modal  energy 
derivatives  from  Eqs.  (31)-(34).  The  first  and  second  electronic  derivatives  were  taken  using  4- 
point  fourth-order  Richardson  Extrapolation11  using  T=.01K  for  h.  They  were  again  taken 
assuming  a  frozen  principal  quantum  number  cutoff  through  the  derivative  steps.  The  absolute 
specific  enthalpy,  hj,  for  each  species  is  calculated  using  Eq.(38)  ,  and  the  energies  previously 
calculated.  The  composition  derivatives  (dCj/3T)p  are  evaluated  numerically  using  Eq.(36).  For 
each  of  the  four  terms,  a  call  to  the  equilibrium  composition  subroutine  is  made  at  the 
appropriate  temperature  (T+AT,  T+ AT/2,  T-AT/2,  T-AT)  keeping  pressure  constant.  The  step 
length  for  the  temperature  is  AT  =  .OIK. 

The  specific  heat  at  constant  volume  is  defined  as 


From  Eq.(35)  we  have  e  =  Ec^,  putting  this  into  the  above 
we  get 
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Unfortunately,  this  cannot  be  easily  evaluated  for  the  following  reason.  The  state  variables  the 
appear  in  the  specific  energy  and  especially  in  the  calculation  of  the  equilibrium  composition  are 
pressure  and  temperature.  A  partial  derivative  can  be  taken  easily  with  respect  to  temperature  or 
pressure  because  we  have  the  ability  to  directly  hold  the  other  constant.  But  here,  we  are 
required  to  take  derivatives  with  respect  to  temperature,  holding  volume  v  (specific  volume) 
constant.  Since  pressure  and  temperature  are  our  chosen  state  variables,  there  is  no  direct  way 
to  do  this.  Authors  such  as  Drellishack4,5,  Peng12,  and  Penski13  use  a  complex  method  involving 
partial  derivatives  of  compressibility,  energy,  and  enthalpy  in  Jacobian-like  form.  But  since  in 
this  research  the  gas  was  treated  as  perfect,  or  at  most  weakly  imperfect,  an  easier  and  more 
fundamental  method  of  obtaining  cw  is  found. 

The  general  relation14 
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is  applicable  to  all  gases,  weather  calorically  perfect,  thermally  perfect,  ideal  or  non-ideal,  and 
even  chemically  reacting.  We  want  to  transform  (de/dv)T  in  to  something  easier  to  calculate.  If 
e=e(v,T)  and  e=e(P,T)  then  it  can  be  shown  that 


Cv  = 


CB~  P 


(46) 


Which,  although  looks  complicated  is  actually  a  much  easier  method  given  the  present  situation 
of  T  and  P  as  state  variables.  cp  is  known,  as  are  the  e^s,  MMG,  and  specific  volume  (from  v 
=  1/p  =  Rugc*T/P*MMG)  because  the  composition  is  already  known  and  the  internal  energies 
have  already  been  calculated  from  the  previous  section.  The  composition  derivative  (0C;/0P)T  is 
evaluated  by  a  4-point  fourth-order  Richardson  Extrapolation  with  AP=.01  N/m2  =  h.  The 
specific  volume  derivative  (0Vj/0T)P  is  evaluated  using  the  same  technique  by  evaluating  the 
molecular  mass  at  the  four  Richardson  step  points  variant  on  T  by 
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with  AT  =  .OIK.  The  derivative  (dMMG/dv)x  requires  the  use  of  an  iterative  technique.  First 
Av  is  calculated  at  each  T  and  P  as  Av  =  .01*v!.  Then  at  each  Richardson  point,  a  vreq  is 
calculated  and  a  variable  Piter  is  guessed  as 

piter  =  pjl_  (48) 


The  equilibrium  composition  subroutine  is  called  using  P=Plter  and  T=T.  A  new  value  of  v, 
called  v“lc  is  obtained  by 
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If  the  absolute  value  of  the  difference  v^-v1^  is  not  less  than  10'9  and  the  maximum  number  of 
iterations  has  not  been  reached  (in  this  case  150),  then  the  pressure  is  updated  to 

J.T  j.»r  ,  RmCT  (50) 

vcalc  v calc  MMG 


and  the  iteration  process  continues.  Convergence  to  required  accuracy  is  usually  obtained  to  from 
within  10-40  iterations.  There  are  some  convergence  problems  at  low  pressures  (P=.001  atm) 
and  high  temperatures  (T> 25,000  K),  but  that  was  rare.  Once  convergence  is  achieved  the 
calculated  molecular  mass  was  set  to  the  appropriate  numerical  derivative  term.  The  specific  heat 
c„  was  then  calculated  using  Eq.  (46). 

The  ratio  of  specific  heats  (y)  is  then  calculated  by 
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which  is  also  its  definition. 


CALCULATION  OF  THE  TRANSPORT  COEFFICIENTS 


In  high-temperature  flows,  the  onset  of  excited  internal  energy  modes,  and  especially  chemical 
reactions  such  as  dissociation  and  ionization,  may  have  a  profound  effect  on  the  way  a  gas 
transports  its  mass,  momentum,  and  energy.  The  microscopic  molecular  transport  of  mass, 
momentum,  and  energy,  arising  from  macroscopic  nonequilibrium  in  composition,  flow  velocity, 
and  temperature,  give  rise  to  the  macroscopic  phenomena  of  diffusion,  viscosity,  and  heat 
conduction,  respectively3.  In  this  research,  the  transport  coefficients  of  viscosity  (n)  and  thermal 


conductivity  (k)  are  calculated  using  the  most  recent  methods  and  data  available.  The  viscosity 
and  thermal  conductivity  are  calculated  including  the  possibility  of  molecular  dissociation,  first 
atomic  ionization,  and  second  atomic  ionization.  The  internal  energy  modes  of  translation, 
rotation,  vibration,  and  electronic  excitation  are  included  when  possible  and  applicable. 
Knowledge  of  the  unique  equilibrium  chemical  composition  and  of  the  equilibrium 
thermodynamic  properties  are  required  for  the  calculation  of  the  transport  coefficients  and  are 
obtained  from  methods  described  in  the  previous  sections. 

A  diatomic  base  gas,  in  this  case  N2,  with  the  inclusion  of  possible  dissociation,  molecular  first 
ionization,  first  atomic  ionization,  and  second  atomic  ionization  gives  six  possible  chemical 
species,  namely  N2,  N2+,  N,  N+,  N++,  and  e'.  These  species  give  rise  to  twenty-one  possible 
independent  binary  interactions,  which  in  this  case  are  illustrated  as  x's  and  o's  in  table  1.  Due 
to  the  lack  of  available  collision  cross  section  data,  and  due  to  the  complexity  of  a  large 
multicomponent  gas  transport  calculation,  some  binary  interactions  and  chemical  species 
exclusions  are  made.  All  interactions  involving  the  molecular  ion  (N2+)  are  ignored.  This  reduces 
the  system  to  a  five  component  transport  model.  The  exclusion  of  N2+  is  based  on  the  fact  that 
it  is  at  all  times  a  minor  species  and  has  negligible  effects  on  the  transport  properties.  The  low 
concentration  is  evident  from  the  previous  equilibrium  composition  calculations.  Peng12  and 
Penski13  made  the  same  assumption.  Capitelli  and  Devoto15  made  a  similar  assertion  but  not  the 
same  assumption,  and  instead  used  crude  estimated  instead  of  absolute  exclusions.  Also,  the 
interactions  involving  N2-N+,  N2-N++,  and  N-N++  were  ignored  because  of  the  low  probability 
that  any  of  these  species  will  co-exist  to  any  significant  degree,  as  evident  from  the  composition 
calculations.  Similar  assumptions  are  made  in  the  works  of  Peng,13  Penski,14  and  Capitelli.15 
These  assumptions  now  leave  us  with  the  twelve  possible  binary  interactions 
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Table  1.  Nitrogen  interaction  model. 


for  the  transport  calculations,  as  indicated  by  an  o  in  table  1.  Although  the  e'-N2  interaction  is 
a  minor  interaction,  it  is  included  because  data  was  available. 

a.  Calculation  of  the  Viscosity 

The  exact  kinetic  representation  for  the  viscosity  of  a  multicomponent  gas  is  given  by 
Hirshfelder16  as 
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In  the  above  formulas  pu  is  the  viscosity  of  a  pure  substance,  given  by 

JmTf 
..  -  /"■  v _ * _ 


=  C 


_2  n(2ar 

a  a  Qjj 


and  /ijj  is  the  binary  gas  viscosity,  given  by 
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Also,  in  the  above  formulas,  Xj's  are  the  mole  fractions  given  by  Eq.(26),  MM;  are  the 
molecular  masses,  and  Ay*  is  the  dimensionless  ratio  of  collision  integrals,  expressed  as 


The  Qy0  s)”s  are  the  dimensionless  collision  integrals 


where  (l,s)  =  (1,1)  corresponds  to  the  diffusion  integral  and  (l,s)  =  (2,2)  corresponds  to  the 
viscosity  and  thermal  conductivity  integral.  They  have  the  physical  significance  of  being  the 
deviation  of  any  particular  model  from  an  idealized  hard  sphere  model.  The  ay's  are  the  collision 
diameter.  And  C  is  a  constant  equal  to  8.4422xl025  when  working  in  MKS  units.  The  only 
difference  between  the  averaged  collision  (transport)  cross  section  and  the  o2Q(l  s)*  is  a  numerical 
constant  equal  to  x  (  o2Q(l,s)*  =  Q/x  ). 16,17 

Calculating  the  viscosity  using  Eq.(52)  is  still  rather  complicated  and  has  proven  to  be 
unnecessary.18  Instead,  Eq.(52)  can  be  replaced  by  its  equivalent  infinite  series  representation, 
which  is19 
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Armaly  and  Sutton19  have  shown  that  the  diagonal  elements  of  matrix  Hy  are  much  larger  than 
the  off  diagonal  elements,  so  much  so,  that  the  off  diagonal  elements  of  Hy  can  be  neglected. 
This  leaves  the  much  simpler  expression  for  the  multicomponent  viscosity  as 
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This  is  modified  by  canceling  an  Xb  such  that 
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Although  the  previous  simplification  seems  trivial,  it  as  actually  quit  computationally  valuable. 
Canceling  out  the  X;  from  each  term  leaves  Hu”s  second  term,  the  summation,  not  dependant 
on  Xj.  This  helps  reduce  the  number  of  0/0  composition  possibilities  in  a  general  purpose 


computational  calculation. 

The  multicomponent  viscosity  of  the  gas  is  then  calculated  using  Eqs.  (61)  and  (62).  The 
composition  is  assumed  to  be  known.  The  data  for  the  collision  integrals  is  sparse,  but  every 
attempt  was  made  to  find  the  most  recent  o2^1,1^  and  o2 Qa,2)*  for  every  reaction.  When  one  or 
the  other  was  unavailable  approximations  had  to  be  made.  Armally  and  Sutton19,20,  Flori  and 
Biolsi21,  and  Freeman  and  Oliver,22  have  all  shown  that  the  ratio  Ay*  is  only  a  rather  weak 
function  of  temperature.  Its  magnitude  is  essentially  governed  by  the  general  type  of  interacting 
particles.  Armally19  gives  the  general  guidelines  that  Ay ’=1.25  for  all  interactions  with  the 
exception  of  atom  interactions  with  its  own  ion.  In  that  case  the  following  values  are 
recommended:  Ay’=0.18  for  H-H+,  Ay*=0.025  for  He-He+,  Ay’=1.10  for  C-C+,  and  for  all 
other  homonuclear  atom-atomic  ion  interactions  the  value  of  Ay’  for  C-C+  should  be  used. 
Freeman  and  Oliver22  have  also  shown  that  although  Ay*  is  a  relatively  weak  function  of 
temperature,  its  temperature  dependance  increases  the  greater  the  pressure.  Thus,  because  of  the 
pressure  range  used  on  this  research,  the  assumption  of  constant  Ay*  was  only  used  when 
absolutely  necessary,  in  this  case  to  get  o2^2,21*  for  the  e'-N  interaction  and  o2^2,21’  for  the  e'-N2 
interaction.  An  interesting  point  is  that  if  Ay’  is  set  equal  to  5/3  in  all  expressions,  then  what's 
left  is  the  mixture  rule  proposed  by  Wilke.  Wilke's  viscosity  mixture  rule  has  been  proven  to  fail 
drastically  for  any  significant  degree  of  ionization.13,19,23  This  happens  for  the  apparent  reason  that 
Ay*  is  never  close  to  the  value  5/3  for  ion-ion  ar  atom-ion  interactions.  The  calculation  of  the 
interactions  with  charged  particles  were  made  using  the  screened  coulomb  potential.24  The  debye 
length  is  taken  as  the  characteristic  dimension  for  coulomb  transport  (  (tW^*  =  XD2Qa,s)* )  which 
in  this  research  is  calculated  as25 

1  (63) 
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in  MKS  units.  The  coulomb  cross  sections  were  calculated  only  if  the  degree  of  ionization  was 
greater  than  10'8.  If  not,  the  cross  sections  were  set  equal  to  a  constant  10'20m2.  This  reduced  the 
effect  of  having  large,  physically  incorrect  values  of  A.D  at  low  degrees  of  ionization.  Although 
calculations  were  done  at  rather  extreme  temperature  and  pressure  conditions,  quantum  effects 
on  the  cross  sections  were  ignored.24,26  The  interpolated  values  of  the  collision  integrals  at  all 
required  points  are  obtained  using  a  cubic-spline  interpolation  technique.  The  calculations  for 
both  the  viscosity  and  the  collision  integral  interpolation  procedure  are  done  in  Fortran. 


b.  Calculation  of  the  Thermal  Conductivity 

The  thermal  conductivity  (k)  for  the  reacting  mixture  is  the  sum  of  three  components,  which 
are  translational,  internal,  and  reactive.  The  basic  translational  relation  for  the  thermal 
conductivity  is  the  familiar  Euken's  relation 
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Anderson2  and  Vincentti  &  Kruger3  give  a  modified  Euken  relation  to  take  into  account  internal 


energy  modes  by  stating  that,  in  general, 
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Also,  for  a  chemically  reacting  gas  ,  there  is  additional  energy  transport  due  to  diffusion,  thus 
giving  as  additional  component  to  the  thermal  conductivity,  called  the  reactive  thermal 
conductivity  (kr).  Since  the  gas  under  consideration  is  homonuclear  with  limited  reaction 
simultaneity,  the  reactive  component  of  the  thermal  conductivity  is  expressed  as12,27 
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where  Dy  is  the  binary  diffusion  coefficient  given  by12 
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Thus  the  governing  equation  for  the  thermal  conductivity  is  now 
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where  kr  is  given  by  Eq.(68). 

For  a  practical  evaluation  of  k,  Eq.(70)  needs  to  be  modified  further.  Recall  from  Section  2 
that  cv  is  a  rather  tricky  quantity  to  get.  The  calculation  of  cv  is  not  accomplished  through  modal 
contributions  (cvttans,  cvrol,  etc.),  but  instead  calculated  as  a  whole  from  chemical  thermodynamic 


relations  and  knowledge  of  cp  using  Eq.(46).  Thus,  there  is  no  general  was  of  obtaining  the 
specific  modal  contributions  and  Eq.  (70)  must  be  modified  as  to  use  the  total  cv.  From  Eq.(64), 
we  have 
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the  total  cv  is 
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now,  re-arranging  the  above  equation  and  putting  this  into  Eq.(70)  we  have 
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Upon  simplification  of  Eq.(73)  and  using  Eq.(64)  again  we  are  left  with 
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where  cv  is  now  the  total  gas  cv.  The  middle  quantity  is  an  equivalent  internal  thermal 
conductivity.  And  kr  is  given  by  Eq.(68). 

The  translational  component  of  the  thermal  conductivity,  k1™5,  appearing  in  Eq.(74)  is  obtained 
from  kinetic  theory.  The  translational  component  of  the  multicomponent  thermal  conductivity  can 
be  expressed  as20,16,21 


£  trans  _ 


h  I*, 

I 

_ jo 


(75) 


where 


L.  =  -4- 


7  trans 

ku 


-  E 


2XtXk 


*- 1  k! 

k*i 


trans 

ik 


(MM, 


*  l 


15  2 

—MM,  + 


-MM;  +  4  MM. MM ^ 


(76) 


and 


2 X. Xj MM. MMj ( 10  -  4 A* ) 
[MMt  +  MMfk^Atj 


where  Xj  is  the  mole  fraction  of  the  i*  component,  Ay*  is  the  dimensionless  ration  of  collision 
integrals  given  by  Eq.(57),  k 7”“  is  the  thermal  conductivity  of  a  pure  gas,  and  ky*”"5  is  the 
thermal  conductivity  of  a  binary  gas.  The  pure  gas  thermal  conductivity  is  given  as 
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and  the  binary  gas  thermal  conductivity  ,ky,  is  given  by 
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where  C  is  a  constant  equal  to  2.6335xl0'23  on  MKS  units.  The  Sly(2,2)*  are  the  dimensionless 
viscosity  and  thermal  conductivity  integrals.  And  ay  is  the  collision  diameter.  The  usage,  theory, 
and  equivalence  to  the  familiar  cross  sections  were  explained  in  the  previous  section. 

As  with  the  viscosity  calculation  (refer  so  Sect.  2)  the  expression  for  ktt  given  by  Eq.(75)  is 
simplified  by  first  expanding  it  in  its  series  representation  as20 
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where  is  has  been  shown  that  the  off  diagonal  elements  Ly  are  much  smaller  than  the  diagonal 
elements.  Thus  a  very  good  approximation  for  the  multicomponent  translational  thermal 
conductivity  is 
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and,  upon  reducing  Eq.(81)  for  the  same  reason  as  done  with  viscosity,  we  are  left  with 
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The  multicomponent  total  thermal  conductivity  is  then  calculated  using  Eq.s(),  (),  and  ().  The 
viscosity  (/*),  the  specific  heat  at  constant  volume  (cv),  the  equilibrium  composition,  and  the 
collision  integral  data  are  assumed  known  (calculated  previously). 


RESULTS  AND  DISCUSSION 

The  electronic  partition  function  and  the  cutoff  principal  quantum  numbers  for  atomic  nitrogen 
are  shown  in  Figs.  1  and  2,  respectively.  Here  we  see  that  due  to  the  cutoff  criterion  of  the 
electronic  partition  function,  a  function  which  is  dependant  on  both  density  and  degree  of 
ionization,  the  electronic  partition  function  is  now  a  function  of  two  state  variables,  namely 
pressure  and  temperature.  The  cutoff  principal  quantum  number  and  the  partition  function  both 
increase  with  increasing  temperature  and  decrease  with  increasing  pressure.  This  general  behavior 
is  present  for  all  other  species,  which  due  to  brevity  are  not  shown.  The  degree  of  ionization, 
shown  in  Fig.  3  illustrates  that  increasing  pressure  retards  ionization.  This  trend  is  also  observed 
for  the  degree  of  dissociation  pictured  in  Fig.  4.  The  total  mixture  thermal  conductivity  is 
displayed  in  Fig.  5.  We  see  that  highly  non-monotonic  variation  occur  in  reactive  regions  where 
the  reactive  component  of  the  thermal  conductivity  dominates.  Beyond  the  dissociation  region, 
although  there  is  some  contribution  due  to  ionization  reactions,  the  translational  component 
dominates.  The  mixture  viscosity,  pictured  in  Fig.  6,  shows  that  the  viscosity  increases  with 
temperature,  independent  of  pressure,  until  reaction  occur  in  the  gas.  The  changes  that  arise  in 
the  presence  of  chemical  reactions  are  due  to  changes  in  the  mixture  molecular  mass  or  in  the 
average  collision  cross  sections.  The  large  increases  that  result  in  the  maxima  are  due  to  the 
dissociation  of  N2.  The  decrease  that  follow  these  maxima  are  due  to  the  increasing  amount  of 
ionization  in  the  mixture.  This  is  caused  by  the  charged  particle  interactions  which  have  large 
coulomb  cross  sections.  The  specific  heats,  given  in  Figs.  7  and  8,  are  strong  functions  of 
temperature,  especially  at  low  pressures  where  the  composition  changes  generally  happen  more 
rapidly  and  the  dc/dT  terms  dominate.  The  ratio  of  specific  heats  is  shown  in  Fig.  9.  It  starts 
at  low  temperatures  at  the  expected  value  of  1.4  (  for  a  diatomic  gas  with  translation  and  rotation 
excited ).  As  the  temperature  increases  we  observe  a  non-pressure-sensitive  rise  due  to  increasing 
vibration  excitation.  As  chemical  reaction  occur,  large  non-monotonic  variations  occur  for  the 
same  reasons  explained  earlier  for  the  specific  heats.  Figures  10  and  1 1  show  that  the  enthalpy 
and  internal  energy  both  increase  with  temperature  and  increase  and  pressure  decreases  in 
reactive  regions  primarily  due  to  the  changes  in  molecular  mass.  Figure  12  clearly  indicates  that 


the  Prandtl  number  is  nowhere  near  constant,  especially  at  high  temperatures,  as  some  frequently 
assume  ( with  as  assumed  constant  value  of  0.75 ).  There  is  a  massive  decrease  with  dissociation 
that  can  cause  it  to  vary  by  a  factor  of  almost  ten. 
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Figure  2.  Cutoff  principal  quantum  numbers  for  N 
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Figure  4.  Degree  of  dissociation  vs.  temperature. 
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Figure  5.  Total  mixture  thermal  conductivity. 
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Figure  7.  Specific  heat  at  constant  volume. 
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Figure  9.  Ratio  of  specific  heats. 
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Figure  10.  Mixture  enthalpy  vs.  temperature. 
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Figure  1 1 .  Mixture  sensible  internal  energy  vs.  temperature. 
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Abstract 


Modeling  of  thermal  diffusion  in  problems  where  the  thermal  conductivity  is  a  severely 
non-monotonic  function  of  temperature  is  studied  numerically.  The  thermal  conductivity  of  a 
partially-ionized,  molecular  gas,  including  transport  due  to  variations  in  chemical  state,  exhibits 
non-monotonic  variations  by  factors  of  three  to  five  in  the  temperature  range  from  1000  -  20,000 
K.  This  range  corresponds  to  conditions  near  surfaces  in  both  internal  and  external  high-speed 
flows,  and  in  various  devices  for  material  processing.  A  detailed  study  of  the  effects  of  numerical 
gridding  on  the  accuracy  of  the  solution  of  the  thermal  diffusion  equation  is  performed  with 
MACH2,  a  2-1/2  dimensional,  time  dependent,  arbitrary  Lagrangian-Eulerian  (ALE)  adaptive- 
grid  magnetohydrodynamic  (MHD)  code.  Data  from  an  advanced  reactive  thermal-conductivity 
model  for  nitrogen  is  provided  in  tabular  form.  An  adaptive  mesh  based  on  variations  in  the 
thermal  conductivity,  rather  than  the  temperature  or  temperature  gradient,  allows  increased 
accuracy  of  calculations  with  a  smaller  number  of  computational  cells.  Comparisons  with  uniform 
gridding  indicate  that,  when  such  an  adaptive  mesh  is  used  with  the  severely  non-monotonic 
transport  coefficient  provided  by  the  advanced  thermal-conductivity  model,  differences  in  thermal 
flux  can  exceed  100%,  especially  at  early  times. 
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Introduction 


In  several  categories  of  aerospace  and  industrial  systems,  partially-dissociated,  partially- 
ionized  gases  are  in  contact  with  solid  surfaces.  These  systems  include  hypersonic  flight  vehicles, 
plasma  thrusters,  and  arc  furnaces.  The  equation-of-state  of  such  a  gas  is  quite  complex  and 
results  in  transport  coefficients  for  thermal  conduction  and  viscous  flow  that  can  be  severely  non¬ 
monotonic  functions  of  temperature  and  pressure.  Solution  of  problems  of  considerable 
engineering  significance,  such  as  heat  transfer  to  ablating  surfaces,  depend  on  accurately 
modeling  transport  across  regions  that  span  a  temperature  range  from  1,000  -  20,000  K.  Over 
this  range,  calculations  for  even  relatively  simple  molecular  gases,  such  as  nitrogen,  indicate  non¬ 
monotonic  variations  of  thermal  conductivity  by  factors  of  three  to  five.  Numerical  modeling  is 
essential  in  order  to  analyze  heat  transfer  under  these  circumstances.  The  question  of  accuracy 
and  efficiency  in  performing  such  modeling  is  the  subject  of  this  paper.  Two  points  of  interest 
are  examined:  1)  the  relation  of  the  calculational  grid  to  the  physical  properties  of  the  problem, 
and  2)  the  use  of  adaptive  gridding  based  on  the  thermal  conductivity,  rather  than  the 
temperature  or  temperature  gradient,  in  order  to  capture  the  effects  of  variation  in  transport 
properties  with  temperature  and  density. 


Example  of  Thermal  Conductivity  Variation 

The  equilibrium  chemical  composition,  thermodynamic,  and  transport  properties  for 
nitrogen  have  been  calculated1.  Six  possible  species  are  included:  N2,  N2+,  N,  N+,  N++,  and  e\ 
The  calculations  are  performed  with  the  inclusion  of  molecular  dissociation  and  vibrational 
excitation,  single  molecular  ionization,  both  single  and  double  atomic  ionization.  Single 
electronic  excitation  is  included  with  a  variable  electronic  partition-function  cutoff  for  atomic 
species  (based  on  distance  between  neighboring  perturbations  to  the  potential  function  of  an 
isolated  atom  or  ion).  The  thermodynamic  state  ranges  for  the  calculations  are:  300  K  <  T  < 
30,000  K  and  10"3  atm  <  p  <  103  atm. 

Figure  1  displays  that  the  thermal  conductivity  of  the  mixture  can  exhibit  highly  non¬ 
monotonic  variations,  especially  in  state  regions  where  chemical  reactions  are  dominant.  For  use 
with  MACH2,  a  2-1/2  dimensional  ALE  MHD  code2,  the  chemical  thermodynamic  and 
transport  data  is  provided  in  the  format  of  SESAME  tables3. 


Numerical  Modeling  of  Unsteady  Heat  Transfer 

It  is  necessary  to  choose  the  time  step  for  numerical  integration  according  to  proper 
characteristic  times  based  on  both  the  temporal  behavior  of  the  boundary  conditions  and  the  cell 
size  needed  to  capture  the  spatial  dependence  of  the  solution.  In  general,  with  the  thermal 
diffusivity  defined  as 
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k 


a) 


K  = 

PCv 

the  characteristic  length  for  unsteady  thermal  diffusion  scales  as 


(2) 


where  tc  is  the  characteristic  time  scale  for  the  boundary  condition,  such  as  the  rise  time  of  a 
sinusoidal  pulse  or  the  elapsed  time  from  the  impulsive  application  of  a  constant  temperature  or 
heat  flux. 

To  capture  the  distribution  of  temperature  within  the  medium,  the  cell  size  must,  of 
course,be  much  smaller  than  this  characteristic  length  scale: 

=  fxxc  (3) 


where  fx  <  <  1. 

The  computational  time  step,  at  least  for  explicit  calculations,  must  be  much  less  that  the 
characteristic  time  for  diffusion  across  the  cell-size: 

6f  =  (4) 

*  K 


where  ft  <  <  1.  Thus,  the  time  step  is  related  to  the  characteristic  time  of  the  problem  by  two 
factors,  fx  and  ft,  such  that: 


8t 


(5) 


Note  that  for  impulsive  application  of  boundary  conditions,  for  which  tc  equals  the  elapsed  time, 
numerical  calculation  will  always  be  inaccurate  at  very  early  times.  Also  the  number  of  cycles 
required  to  complete  computation  of  a  problem  of  duration  tc  is 


N  = 


(6) 


For  ft=fx=  0.1,  N  =  1000  cycles.  This  is  overcome  by  the  use  of  an  implicit  scheme. 
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In  computations  presented  here,  in  which  the  cell  size  is  changing  throughout  the  computation, 
the  maximum  time  step  allowed  is  chosen  to  be  equal  to  the  characteristic  time  step  to  minimize 
the  influence  of  errors  that  might  otherwise  result  from  improperly  chosen  time  steps  for  the 
smallest  cell,  with  f,=  l. 


One-Dimensional  Thermal  Flux  on  a  Uniform  Grid 

Computations  were  performed  for  three  values  of  evenly-spaced  cells  (n  =  64,  32,  and 
16),  respectively,  Figure  A  displays  the  temperature  profile  for  unsteady,  one-dimensional,  heat 
transfer  between  two  surfaces  impulsively  imposed  across  a  medium  initially  at  a  uniform 
temperature  (equal  to  the  average  value  of  the  two  surfaces)  with  a  constant  thermal  conductivity. 
Figure  B  compares  the  error  in  calculations  of  temperature  as  a  function  of  time  between  the 
surfaces,  at  several  times,  in  terms  of  the  differences  between  local  values  of  temperature  (for 
n  =  64,  32,  and  16,  respectively)  and  the  local  value  obtained  analytically4.  This  error  is 
normalized  by  the  mean  value  and  given  as  a  cell  averaged  quantity  across  the  grid.  This  is 
calculated  as: 


no. cells 


xlOO  <7) 


Accuracy  at  early  times  improves  as  the  cell  size  becomes  smaller  compared  to  the  characteristic 
length  for  diffusion.  Smaller  cell  size,  however,  within  the  strict  requirement  of  Eqn.  4,  requires 
a  larger  number  of  computational  cycles.  This  requirement  is  certainly  ameliorated  by  the  use 
of  implicit  solution  of  the  difference  equations,  but  at  the  potential  loss  of  accuracy  even  within 
stable,  converged  solutions.  The  use  of  nonuniform  cell  size  permits  a  reduction  in  the  total 
number  of  calculations  and  memory  requirements,  if  the  reduced  number  of  cells  can  be 
concentrated  in  regions  where  important  variations  are  occurring. 


Adaptive  Gridding 

The  benefit  of  using  a  calculational  grid  with  variable  resolution  is  well  established. 
Often,  the  grid  will  be  determined  by  the  shape  of  flow  boundaries  and  by  the  distribution  of 
flow  functions  such  as  mass  density.  Figures  4  and  5  display  the  percent  differences  in  flux  of 
the  solutions  for  a  unsteady  heat  transfer  problem  using  the  tabular  conductivity  with  a  grid  that 
is  concentrated  in  regions  of  highest  temperature  gradient  to  that  of  an  ordinary  fixed  grid.  This 
is  done  for  16  and  32  cells,  respectively. 

Here  we  see  that  the  largest  differences  happen  at  earlier  times,  especially  near  the  domain  mid¬ 
range,  and  decrease  as  time  continues.  Larger  differences  appear  in  the  32  cell  case.  The 
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magnitudes  being  approximately  double  that  of  the  16  cell  simulation. 


Grid  Based  on  Thermal  Conductivity 


The  code  has  the  capability  to  adapt  the  grid  according  to  the  magnitude  of  the  thermal 
conductivity  and,  separately  or  in  combination,  to  gradients  in  the  thermal  conductivity.  The 
decision  to  have  the  ability  to  adapt  the  grid  on  either  magnitudes  or  gradients  is  based  on  the 
assumption  that  both  may  be  important  depending  on  the  physical  situation  or  the  transport 
model. 


The  specific  amount  of  adaptation  provided  be  each  control  is  determined  by  the  following 
procedure.  This  is  accomplished  by  continuously  adjusting  the  grid  to  the  squares  of  the 
logarithmic  derivatives,  accumulated  of  course,  throughout  the  spatial  and  temporal  coordinates. 
The  exact  sensitivity  of  the  calculational  grid  is  defined  by  user  defined  input  parameters  which 
modify  the  code  internal  adaptive  weight  function  accordingly. 

The  ability  to  adapt  according  to  the  thermal  conductivity  in  this  manner  has  four  major 
features.  Enough  cells  are  placed  in  regions  where  large  amounts  of  thermal  diffusion  are 
important.  The  true  maxima/ minima  of  conductivity  are  captured  in  descretized  space  when 
dealing  with  highly  non-monotonic  transport  model,  thus  minimizing  the  risk  of  missing 
important  values  by  approximations  in  using  the  tabular  data.  The  thermal  conductivity  is  a 
function  of  both  pressure  and  temperature,  so  the  grid  directly  adapts  on  one  function  instead  of 
two,  thus  minimizing  computational  work.  Finally,  the  exact  correlation  of  the  transport  model 
with  flow  structure  can  be  delineated. 

To  test  the  adaptive-gridding  scheme,  the  unsteady,  two-wall  heat  transfer  problem  is 
again  considered.  As  pictured  in  Fig  2,  the  left  wall  is  at  1.0  ev,  the  right  at  0.5  Ev,  and  the  gas 
is  initially  at  a  uniform  temperature  of  0.75  Ev.  Before  using  the  more  complex  tabular  model 
(for  Fig.  1),  it  is  useful  to  employ  a  generic  ersatz,  analytical  model  for  a  severely  non¬ 
monotonic  variation  in  thermal  conductivity.  A  sinusoidal  form  is  chosen  for  this  task.  To  test 
the  adaptive  grid  algorithm,  we  impose  a  sinusoidally- varying  thermal  diffusivity  on  the  interval 
[0,1]: 


K(x)  =  1+  bsm.(2mtx)  ] 


(8) 


For  Kfl  =  0.25  m2  /  s  ,  b  =  0.5,  and  n  =  4,  the  resulting  conductivity  profile  is  illustrated  in 
Figure  3. 
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Twelve  different  runs  were  performed  utilizing  the  sinusoidal  thermal-conductivity  model 
given  by  Eqn.  6,  under  four  different  adaptivity  conditions.  Three  runs  using  a  plain  non- 
adaptive  grid  were  performed  as  baseline  computations  using  64,  32,  ad  16  cells,  respectively. 
The  results  of  the  64  cell  run  are  shown  in  Fig.  8.  The  results  of  the  32  and  16  cell  runs 
exhibited  similar  behavior.  The  64  cell  grid  simulation  using  the  adaptive  grid  based  on  thermal 
conductivity  magnitude  is  shown  in  Fig  9.  The  grid  density  for  this  case  is  shown  in  Fig.  10. 

Figures  11  and  12  show  the  percent  difference  for  the  sinusoidal  conductivity  model 
between  the  uniform  grid  and  the  grid  adapting  to  the  magnitude  for  the  16  and  32  cell  cases, 
respectively.  The  difference  being  defined  as: 


E(xa,ta) 


Qg(XM  -  Qf^gSg) 


X  100 


(9) 


where  Qf(xa,ta)  is  the  interpolated  value  of  flux  from  the  fixed  grid  simulation  at  each  adaptive 
grid  point.  (This  interpolation  procedure  is  performed  using  a  multi-dimensional  cubic  spline 
interpolation  ).  The  differences  indicated  are  very  small  except  in  the  mid-range  where  a  rather 
large  variation  occurs.  For  the  32  and  16  cell  cases,  much  larger  errors  exist  near  the  middle 
reaching  consistent  maxima  of  approximately  78%  for  the  32  cell  case  and  66%  for  the  16  cell 
case.  Note  the  small  deviation  even  at  late  times  for  both  of  these  cases. 

Calculations  with  the  true  non-monotonic  conductivity  model  of  Fig.  1,  under  the  same 
general  parameters  as  the  simulations  with  the  idealized,  sinusoidal  form,  yield  similar  results. 
Figure  13  shows  the  flux  distribution  with  no  adaptivity.  Figure  14  shows  the  resultant  flux 
distribution  with  full  adaptivity.  The  percent  difference  between  the  fixed  grid  and  fully  adaptive 
grid  for  64  cells  is  shown  in  Fig.  15.  We  see  that  the  differences  present  in  the  tabulated  runs 
are  generally  much  greater  that  present  in  the  sinusoidal  model.  Figure  15  shows  a  rather  large 
maxima  in  the  mid-range  with  a  maximum  of  approximately  28%.  Note  that  the  left  wall  maxima 
is  still  present  in  this  case  along  with  the  32  and  16  cell  cases  shown  in  Figs.  16  and  17.  Figures 
18  through  20  show  the  grid  densities  for  these  three  cases. 


Concluding  Remarks 

Calculation  of  a  multi-component,  chemically-reacting  transport  model  for  nitrogen 
revealed  that  in  certain  state  regimes  the  thermal  conductivity  exhibits  highly  non-monotonic 
behavior1.  This  prompted  the  modification  of  the  MACH2  computer  code  to  adapt  the 
computational  grid  to  the  magnitude  and  gradient  of  the  thermal  conductivity  and  the  gradient 
of  the  thermal  conductivity.  As  time  progresses,  such  an  adaptive  grid  increases  the  accuracy  of 
the  solution  by  concentrating  the  grid  in  the  appropriate  regions  to  capture  severe  variations  in 
thermal  conductivity.  This  adaptation  also  helps  to  stabilize  the  numerical  solution  at  boundaries 
with  impulsive  initial  conditions.  The  largest  percentage  differences  between  results  with  the 
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new  gridding  scheme  and  uniform  gridding  occur  at  early  times,  but  errors  persist  even  at  late 
times  in  the  calculations. 

In  addition  to  avoiding  significant  errors  in  calculation  of  heat  flux  (in  excess  of  70%  in 
some  instances),  the  use  of  a  gridding  scheme  that  follows  the  behavior  of  the  thermal 
conductivity  will  be  important  in  attempts  to  understand  the  structure  of  nonuniform  regions 
between  relatively  cold  boundaries  and  high  temperature,  molecular  flows.  Aerospace  problems, 
including  heat  transfer  and  ionization  in  hypersonic  boundary-layers,  and  the  current  distribution 
in  ablation-fed  plasma  thrusters,  will  certainly  depend  on  accurate  modeling  of  transport 
properties,  but  also  require  that  numerical  computations  capture  the  complexity  of  such  modeling. 
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Figure  1.  Total  mixture  thermal  conductivity  for  nitrogen. 


Temperature  Distrubution 

k  =  const.  ,  si  =s2=0 


1.269e-3 
0.1002  s 
0.2005  s 
0.3007  s 
0.401  0  s 
0.5012  s 
0.6002  s 
0.7004  s 
0.8007  s 
0.9001  s 
1 .001  2  s 


0.501 - * — 

0.0  0.1 


0.5 
X  (m) 


Figure  3.  Temperature  distribution.  No  adaptivity.  Constant  conductivity. 
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Figure  4.  Error  between  Mach2  calculations  and  analytic  solution. 


Figure  5.  Percent  difference  in  flux  for  16  cells.  Tabular  conductivity. 
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ABSTRACT 

The  chemical  composition  and  thermodynamic  properties  of  tetrafluoroethylene  (C2F4)  is 
calculated  with  a  two-temperature  LTE  formulation.  Twenty-five  chemical  species  are  included 
in  the  analysis.  The  equilibrium  constants  are  calculated  using  the  most  recent  spectroscopic  data 
available.  Calculations  are  performed  for  pressures  from  0.001  atm  to  1.0  atm  and  for 
temperature  ranges  of  0.05  ev  to  5  ev  for  both  heavy  particle  and  electron  temperatures. 


INTRODUCTION 

Knowledge  of  the  chemical,  thermodynamic,  and  transport  properties  of  a  gas  is  required  in 
almost  any  gasdynamic  analysis.  Accurate  thermochemical  and  transport  properties  become 
particularly  important  in  high-temperature  applications  such  as  the  pulsed  plasma  thruster.  In  this 
paper,  we  shall  concentrate  on  calculating  the  equilibrium  composition  of  a  gas  mixture.  This 
is  the  necessary  first  step  for  determining  the  thermodynamic  and  transport  properties  of  a  gas. 

There  are  three  primary  thermochemical  states  possible  for  a  gas.  A  calorically  perfect  gas  has 
specific  heats  that  are  constant,  and  the  enthalpy  and  internal  energy  are  only  functions  of 
temperature.  A  thermally-perfect  gas,  in  which  variable  vibrational  and  electronic  excitation  are 
taken  into  account,  has  specific  heats,  enthalpy,  and  internal  energy  that  are  all  functions  of 
temperature.  If  the  conditions  are  right  for  chemical  reactions  to  occur,  then  we  can  treat  the  gas 
as  an  equilibrium  chemically-reacting  gas  for  which  properties  area  all  functions  of  temperature 
and  pressure.  Even  this  can  be  generalized  by  stating  that  the  reacting  gas  is  in  local 
thermodynamic  equilibrium  (LTE).  This  means  that  a  local  Boltzmann  distribution  exists  at  each 
point  in  the  flow  at  the  local  temperature.  We  will  extend  this  statement  further  for  the  case  of 
a  two  temperature  LTE  gas  modeled  here.  In  this  paper  we  calculate  the  chemical  composition 
of  gaseous  Teflon. 

Very  little  information  exists  on  the  thermochemical  properties  of  PTFE  and  its  chemical 
constituents,  especially  in  the  plasma  state.  This  is  unfortunate  for  high-temperature 
thermochemical  and  transport  data  is  needed  for  such  devices  as  pulsed  plasma  thrusters1-2, 
ablative  re-entry  heat  shields,  electronic  components  (such  as  capacitors),  plasma  etching 
systems,  and  upper  atmosphere  modeling.  Nonequilibrium  properties  are  particularly  important 
when  studying  electric  discharges  near  Teflon  surfaces  and  other  nonequilibrium  flows. 


POSSIBLE  SPECIES,  REACTIONS,  AND 
EQUILIBRIUM  EQUATIONS 


In  this  paper,  we  calculate  the  chemical  composition  of  tetrafluoroethylene  (C2F4).  The  analysis 
will  include  vibrational  and  electronic  excitation,  dissociation,  first  molecular  ionization,  and  first 
through  fourth  monatomic  ionization.  Throughout  of  the  analysis,  we  shall  assume  a  perfect  gas, 
where  intermolecular  forces  are  non-existent  or  negligible.  This  might  seem  a  strange  assumption 
when  the  gas  is  in  the  plasma  state  due  to  the  presence  of  Coulomb  collisions,  but  it  is  a  widely 
used  and  accepted  approximation.3 

For  a  polyatomic  base  gas,  C2F4  in  our  case,  with  the  possibility  of  undergoing  full 
dissociation,  singular  molecular  ionization,  and  up  to  fourth  monatomic  ionization,  we  first 
assume  there  are  twenty-five  possible  chemical  species,  which  are  C2F4,  C2F2,  CF2,  CF2+,  CF3, 
CF3+,  CF4,  C2,  CF,  CF+,  F2,  F2+,  Cz  (Z=-l,4),  Fz  (Z=-l,4),  and  e*  's.  For  a  gas  containing 
twenty-five  chemical  species,  which  is  composed  of  three  elements  (C,  F,  e  ),  we  are  required 
to  have  twenty-two  (25-3=22)  independent  chemical  reaction  equations  (laws  of  mass  action). 
The  reactions  considered  here  are 


CF2  *±CF  +  F  (1) 

CF  *±C  +  F  (2) 

F2^2F  (3) 

CF2  ^  CF2  +  e~  (4) 

CF  **  CF*  +  e ~  (5) 

Cz  1  **  Cz  +  e~  ]  (6) 

Z=0 

C2F4^2CF2  (7) 

F2  ^  f;  +  e-  (*) 

Fz-i  _  Fz  +  e-  |  (9) 

Z= 0 


C2F4  ^  C2F2  +  f2 

(10) 

C2Fa  -  CF3  +  CF 

(ID 

c2fa  **  CF4  +  c 

(12) 

c2  ^  2  C 

(13) 

In  actuality,  there  are  other  possible  reactions  that  could  yield  the  same  chemical  species.  But, 
for  an  equilibrium  calculation,  the  reactions  chosen  are  arbitrary  as  long  as  they  are  linearly 
independent  and  account  for  all  possible  species. 

Writing  these  reactions  in  terms  of  equilibrium  relations  for  the  partial  pressures,  we  have 

=  n  p,‘  <»4> 


where  the  Kpj  are  the  equilibrium  constants  for  the  reaction  (j)  at  the  equilibrium  temperature 
T,  in  terms  of  the  partial  pressures.  Using  the  appropriate  formulations,  they  may  also  be  put 
in  terms  of  concentrations,  Kc,  or  number  densities,  Kn.  It  is  important  to  note  in  the  above 
equations  that  the  equilibrium  constants  are  written  as  functions  of  temperature  only,  as  most 
authors  point  out.  However,  they  may  be  functions  of  two  or  more  state  variables  depending  on 
whether  such  things  as  thermal  non-equilibrium  assumptions  or  electronic  partition  function 
cutoff  is  taken  into  account.4 

In  addition  to  twenty-two  independent  equations  relating  the  twenty-five  unknown  partial 
pressures,  we  need  three  more  equations  to  solve  for  the  gas  composition.  The  three  chosen  are; 
conservation  of  nuclei,  Dalton’s  Law,  and  charge  neutrality.  The  ideal  thermal  gas  law  for  each 
species  is  written  in  the  form 

P ,  =  »,*?■,  (15) 

where 

p  =  E  Pi  (16) 

For  charge  neutrality,  we  have 
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In  terms  of  partial  pressures,  this  becomes  (  for  ideal  gases  ) 
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For  conservation  of  nuclei,  we  write 
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where  (nC2F4)o  is  the  total  number  of  tetrafluoroethylene  molecules  available  for  dissociation  and 
ionization  (ie.  the  number  of  C2F4  molecules  present  if  the  gas  was  non-reacting  at  some  initially 
low  temperature).  Dividing  Eq.  (20)  by  Eq.  (21)  and  utilizing  Eq.  (19)  gives  us  the  nuclei 
conservation  statement,  where  the  number  densities  are  related  to  the  partial  pressures  by  Eq. 
(15). 


CALCULATION  OF  THE  EQUILIBRIUM 
CONSTANT  -  PARTITION  FUNCTIONS 


To  solve  the  system  of  equations,  we  only  need  values  for  the  equilibrium  constants  which 
may  be  calculated  from  equilibrium  statistical  mechanics.  In  terms  of  partition  functions  Qh  the 
law  of  mass  action  for  a  general  system  is 


=  n  n'.:  - 


(12) 


or  alternatively,  substituting  n;  =  N/V  we  have 
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where  V;  is  the  stoichiometric  mole  number  for  species  (i),  that  is,  the  coefficients  in  the 
balanced  chemical  equation,  Ae0  is  the  reaction  energy  (change  in  zero-point  energy)  and  is 


the  total  partition  function  for  species  (i).  Thus,  for  a  given  reaction  and  thermodynamic  state, 
the  only  unknowns  in  Eq.  (23)  are  the  Q;’s. 

For  a  system  in  thermodynamic  equilibrium,  we  have 


(24) 


which  gives  the  number  of  particles  N/  in  energy  level  ti  with  gj  degenerate  states.  We  define 
the  partition  function,  Q,  as  the  sum  in  the  denominator  of  Eq.  (24) 


Q  -  Ear/1* 


(25) 


which  is,  in  general,  a  function  of  T  and  V.  It  is  typical  to  express  the  total  energy  as  the  sum 
of  translational  and  internal  energies.  Note  that  Eqs.  (23)  and  (24)  contain  only  one  temperature. 
For  the  two  temperature  case  considered  in  this  research  we  make  the  assumption  that  the  heavy- 
particle  gas  composed  of  neutrals  and  ions,  has  a  Maxwellian  distribution  in  velocities  and  a 
Boltzmann  distribution  in  energies  at  a  heavy-particle  temperature  ,  T.  The  electron  gas, 
composed  of  both  free  and  bound  electrons  is  in  equilibrium  with  an  electron  temperature  Te 
defined  by  their  Maxwellian  velocity  distribution.  Note  that,  in  this  analysis,  we  are  ignoring  the 
interaction  between  electronic  and  vibrational  states.  Thus,  we  have  defined  a  two-temperature 
LTE  situation.  Equation  (22)  is  the  standard  statistical  representation  of  the  equilibrium  constant 
for  a  one  temperature  system.  For  the  multi-temperature  system  analyzed  here  we  expect  a 
different  representation.  There  has  been  some  controversy  about  the  correct  form  of  the  multi¬ 
temperature  Saha  equation5,6,7,8.  Derivations  using  kinetic  arguments  yield  the  same  form  as  Eq. 
(22)  but  with  T  replaced  by  Te  in  the  Boltzmann  factor.  This  formulation  has  long  been 
consistent  with  the  general  assumption  that  ionization  is  in  equilibrium  with  the  electron 
temperature  and  just  replacing  T  by  Te  as  an  approximation.  Other  derivations  using  pure 
thermodynamic  arguments  obtain  a  Saha  equation  in  which  both  temperatures  appear.  This 
inconsistency  was  examined  in  detail  by  Van  de  Sanden5,  et  al.  and  was  found  that  these 
representations  were  derived  using  the  standard  form  of  the  generalized  law  of  mass  action,  that 
is 


E  l*.».  *  «  (26) 


which  is  proven  incorrect  for  a  multi-temperature  system  because  it  may  violate  the  second  law 
of  thermodynamics  for  the  system  as  a  whole.  Van  de  Sanden  found  the  correct  form  to  be 
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Which  when  used  in  the  thermodynamic  derivation  results  in  a  Saha  equation  equivalent  to  that 
obtained  by  the  kinetic  one.  Thus  in  this  research  the  equilibrium  constant  for  ionization  reactions 
is  Eq.  (22)  with  Te  in  the  exponential  term. 

For  a  molecule  we  have 
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and  for  an  atom 
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where  e  is  the  sensible  energy,  measured  above  the  zero-point  energy  e0.  Quantum  mechanics 
gives  us  theoretical  values  for  the  quantized  energies  of  a  particle,  at  least  for  the  translational, 
rotational,  and  vibrational  modes.9  Along  with  the  assertion  that  particle  energy  is  simply  the  sum 
of  the  modal  energies,  that  is,  the  internal  energies  are  uncoupled,  which  is  a  consequence  of  the 
more  fundamental  assumption  of  a  separable  Hamiltonian,  the  partition  function  is  expressed  as 
the  product  of  the  modal  partitions  Qj5  where 

<?  =  !!<?,  <30> 


with  j  extending  over  all  modes,  that  is,  translational,  rotational,  vibrational,  and  electronic 
modes.  Armed  with  the  quantized  values  for  the  modal  energies,  and  the  associated  degeneracies 
we  can  calculate  the  modal  partition  functions  which  are  given  here  in  reduced  form,  without 
proof  as10 
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where  a  is  a  factor  which  arises  from  the  symmetry  requirements  of  the  wave  function  in  the 
exchange  of  an  identical  particle.  It  is  equal  to  1  for  heteronuclear  molecules  (ex.  CF),  and  equal 
to  2  for  homonuclear  molecules  (ex.  F2).  This  number  may  take  other  values  depending  on  the 
complexity  of  the  molecular  structure.  The  rotational  and  vibrational  partition  functions  are 
derived  assuming  that  the  molecule  is  a  rigid  rotator  and  harmonic  oscillator,  respectively.  For 
any  general  polyatomic  molecule  of  N  atoms,  if  we  still  assume  a  separable  Hamiltonian  then 
we  can  factor  the  partition  function  as  in  Eq.  (30)  with  the  product  extending  not  only  over  all 
fundamental  modes  but  also  over  all  modal  degrees  of  freedom.11  For  a  polyatomic  molecule 
composed  of  N  atoms,  there  is  a  total  of  3N  geometric  degrees  of  freedom.  Translation  requires 
3  degrees  of  freedom  thus  there  are  a  total  of  3N-3  degrees  of  freedom  available  for  the  internal 
modes.  For  linear  polyatomic  molecules  the  rotation  can  be  completely  specified  by  two  angles, 
thus  taking  up  two  more  degrees  of  freedom.  This  leaves  3N-5  degrees  of  freedom  for  vibration. 
The  rotational  partition  function  if  still  given  be  Eq.  (32).  For  nonlinear  molecules  three  angeles 
are  necessary  thus  there  are  a  total  of  3N-6  geometric  degrees  of  freedom  for  vibration.  The 
rotational  partition  function  in  this  case  is  given  by11 
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For  each  case  the  total  vibrational  partition  function  is  taken  as  the  product  of  each  of  the 
vibrational  degrees  of  freedom  with  the  partition  function  for  each  degree  of  freedom  given  by 
Eq.  (33).  As  stated  earlier,  we  have  assumed  uncoupled  energy  modes.  That  is,  we  assume  that 
the  energy  exchange  between  energy  modes  is  small.  As  a  result  each  energy  mode  which  is 
equilibrium  at  its  temperature,  in  our  case  either  T  or  Te,  is  statistically  independent  from  all 
others.  Thus  for  our  multi-temperature  system  we  can  write 

q  =  n  (35) 


where  T;  is  the  equilibrium  temperature  of  the  mode. 

Note  that  the  Boltzmann  factor,  e‘Bej,  in  the  above  equations  for  the  partition  functions  is 
written  as  derived  in  the  standard  one-temperature  system  case.  That  is,  where  8  =  1/Kt.  It  has 
been  known  for  a  long  time  that  certain  molecules  have  the  tendency  to  have  their  vibrational 
modes  to  be  a  strong  function  of  electron  temperature,  that  is  Tv  =  Te  instead  of  T.  Atmospheric 
gasses,  especially  nitrogen  are  well  known  for  this  property12,13.  This  phenomena  is  thought  to 
be  largely  due  to  the  propensity  of  the  molecule  for  vibrational  resonance.  Since  we  are  dealing 
with  a  gas  mixture  containing  twelve  molecular  species  a  more  formal  examination  of  the 
situation  is  required  to  support  our  assumptions.  We  define  a  parameter  r\  which  quantifies  the 
fractional  energy  exchange  between  thermal  systems,  in  our  case  it  is  the  fractional  energy 
exchange  between  systems  at  the  equilibrium  electron  temperature,  ie.  the  electron  gas,  and  other 
internal  modes.  From  the  results  derived  in  Appendix  1  we  find  that  our  Lagrange  multiplier  B 
becomes 
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At  this  point  in  time,  research  is  being  done  on  how  to  set  the  ti  parameter  effectively. 
Consistent  with  the  assumptions  made  thus  far  we  have  for  the  interaction  of  free  electrons  and 
bound  electrons  rj  =  1  thus  6  =  1/Kte  for  Qe.  And  for  the  interaction  of  free  electrons  and  the 
translational,  rotational,  and  vibrational  modes  q  =  0  thus  8  =  I/K4  for  the  partition  functions 
for  these  modes. 

For  electronic  energy  there  is  no  closed  form  general  expression  for  the  quantized  energy 
levels,  thus  the  electronic  partition  function  must  be  left  as  an  infinite  series  in  the  form 


j=0 


(37) 


Equation  (37)  is  the  true  theoretical  representation  of  the  electronic  partition  function  for  an 
isolated  particle.  In  theory  there  are  an  infinite  number  of  electronic  levels  extending  from  the 
ground  state  energy  (e0  =  0)  to  the  ionization  potential,  which  is  the  amount  of  energy  needed 
to  remove  an  electron  from  its  ground  state  to  infinity  ( ie.  undergo  a  bound-free  transition).  The 
electronic  partition  function  is  a  diverging  series  because  although  the  energy  approaches  a  finite 
limit,  the  degeneracy  increases  as  the  square  of  the  principal  quantum  number,  so  the  series 
diverges.14 

In  actuality,  the  electronic  series  in  not  infinite  because  a  particle  in  the  real  world  is  never 
truly  isolated.  Due  to  various  interparticle  interactions  that  arise  in  any  finite  density  medium, 
the  series  will  actually  terminate  at  some  principal  quantum  number,  ncutoff.  The  evaluation  of 
ncutoff  and  its  associated  effect  of  ionization  potential  lowering  is  the  subject  of  some  controversy 
and  was  explored  in  detail  in  another  work.4  Results  from  that  work  give  the  correct  cutoff 
criterion  as 
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where  2— Zeff*  e2/2*IP.  The  lowered  ionization  potential  is  given  by 
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The  K„'s  are  converted  to  Kp's  using  the  relation 
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The  required  molecular  and  atomic  data,  which  are  too  numerous  to  give  here,  are  taken  from 
the  works  of  Chase,15  Moore,1617  Rosenstock,18  Herzberg,19  Buckely,20  and  Paulino  and  Squires.21 


SOLUTION  OF  THE  EQUILIBRIUM 
EQUATIONS 

The  equations  given  in  the  previous  section  to  calculate  the  equilibrium  composition  in  terms 
of  partial  pressures,  given  T  ,  Te  and  P,  provide  a  closed  set  of  twenty-five  coupled  nonlinear 
algebraic  equations  for  which  there  is  no  analytic  solution  and  numerical  methods  must  be  used. 
The  numerical  solution  of  systems  of  nonlinear  equations  is  universally  very  difficult  and  is  a 
topic  of  current  research.  As  of  yet,  there  are  no  appropriate  numerical  methods  for  solving 
coupled  nonlinear  systems  of  algebraic  equations  from  arbitrary  starting  vectors.22  One  of  the 
most  commonly  used  methods  and  the  one  used  previously  for  nitrogen  research,4  is  the  Newton- 
Raphson  method.  The  Newton-Raphson  method  usually  exhibits  excellent  convergence  qualities 
when  the  starting  vector  (initial  guesses  for  the  roots)  is  near  the  actual  root.  For  a  homonuclear 
diatomic  gas  a  good  initial  guess  can  usually  be  obtained  through  Saha-type  statistical  arguments 
using  weighted  averages  due  to  limited  reactive  simultaneity.  In  this  case  though,  since  we  are 
dealing  with  a  polyatomic  heteronuclear  gas,  it  was  initially  assumed  that  we  could  not  get  close 
enough  to  the  root  for  Newton-Raphson  to  work  properly.  Many  hybrid  techniques  have  been 
proposed  but  most  fail  when  the  Jacobian  becomes  singular,  or  at  stationary  point.23,24  The 
solution  procedure  initially  chosen  for  this  research  is  the  one  proposed  by  Powell.25  Powell's 
technique  exhibits  almost  guaranteed  convergence  even  for  poor  initial  guesses.  It  also  has  the 
ability  to  correctly  handle  stationary  points.  The  price  of  this  behavior  is  that  the  convergence 
is  linear  until  very  close  to  the  root  then  it  converges  quadratically  like  the  standard  Newton- 
Raphson.  Thus,  the  total  number  of  iterations  required  is  greatly  increased  resulting  in  much 
longer  runtimes,  up  to  200  minutes  per  isobar  on  a  workstation  .  Thus,  it  was  decided  that  the 
Newton-Raphson  method  would  be  used  combined  with  a  raster  processing  iteration  procedure 
throughout  the  computational  thermodynamic  space. 

Given  T,  Te,  and  P,  an  initial  equilibrium  composition  is  approximated  using  simple  Saha-type 
arguments  assuming  totally  uncoupled  physical  process  at  low  T  and  Te,  then  partial  pressures 
are  backed  out  using  weighted  approximations  or  the  initial  guess  was  set  equal  to  the  previous 
converged  root  at  the  last  temperature  and  pressure  and  stepped  up  in  temperature  on  each  isobar 
in  small  increments  of  from  0. 1  to  10  K  for  heavy  particle  temperatures  and  1  K  to  100  K  for 
electron  temperatures.  The  completely  closed  system  is  solved  at  each  point  using  the  Newton- 
Raphson  procedure,  which  iterates  until  the  sum  of  the  absolute  values  of  the  corrections  is  less 
than  the  chosen  tolerance,  C,  where 
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where  e  =  10-4  or  until  a  maximum  of  20,000  iterations  are  performed.  Convergence  to  the 
tolerance  is  usually  obtained  quickly.  An  additional  convergence  check  is  performed  after  each 
completed  iteration  so  that  the  total  isobaric  deviation  is  not  allowed  to  exceed  5%.  The  correct 
composition  is  now  known  to  within  desired  accuracy.  The  mixture  molecular  is  computed  using 
the  expression 
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where  X;  is  the  mole  fraction  of  species  i,  obtained  using 
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where  n^  is,  of  course,  just  the  total  number  of  particles. 


CALCULATION  OF  THE  THERMODYNAMIC  PROPERTIES 


For  a  given  temperature  and  pressure,  the  complete  gas  composition  is  calculated  using  the 
techniques  in  the  previous  section.  Often,  however,  it  is  necessary  to  have  the  thermodynamic 
properties  of  the  gas.  In  this  research,  the  calculated  thermodynamic  properties  are  specific 
enthalpy  (h),  specific  internal  energy  (e).  The  specific  heats  (Cp  and  cv),  and  consequently  the 
ratio  of  specific  heats  (y),  were  not  calculated  for  there  is  some  ambiguity  as  to  the  definition 
of  specific  heat  for  a  multi-temperature  system  and  they  are  rarely  used  in  high-temperature 
calculations  anyway. 

On  a  microscopic  scale,  a  particle  (atom  or  molecule)  may  possess  energy  due  to  its 
translational,  rotational,  vibrational,  and  electronic  modes.  This  is  expressed  as 
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where  e  is  the  total  specific  internal  energy.  Note  that  e  is  actually  the  total  specific  sensible 
internal  energy  and  is  measured  above  the  zero  point  (e=e-e0-i  where  eo  i  is  the  summation  of  the 
zero  point  energies  for  the  individual  modes  of  translation,  vibration,  and  electronic.  The 
rotational  mode  has  no  zero  point  energy.  In  terms  of  partition  functions,  the  energy  and 
enthalpy  for  a  one  temperature  system  is  expressed  as10 


(45) 
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with  the  relation 

h  =  e  +  Pv  (47) 


holding  for  any  type  of  gas,  whether  calorically  perfect,  thermally  perfect,  or  chemically 
reacting.  The  expressions  for  the  modal  partition  functions  explained  earlier.  Taking  the 
derivatives  of  the  logarithms  and  using  Eq.  (45)  we  get 
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Again,  for  the  electronic  energy,  the  lack  of  a  closed  form  partition  function  leaves  us  with 
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where  Rj  =  RUGC/MMi  is  the  specific  gas  constant.  The  translational  and  rotational  modes  are 
assumed  to  be  fully  excited.  For  a  multi-temperature  system,  composed  of  n  subsystems  each  in 
equilibrium  at  their  own  temperature  Tn,  as  long  as  the  subsystems  are  statistically  independent 
then  the  fundamental  thermodynamic  relations  may  be  applied  to  each  temperature  subsystem 
separately5,  thus  we  write  Eqs.  (45)  and  (46)  as 
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For  a  mixture  of  N  species,  the  total  mixture  sensible  internal  energy  is 


N 


e  =  E  ci«. 


i  =  l 


where  q  is  the  mass  fraction  of  species  i,  given  by 

mi  P, 
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m  p 


Similarly,  the  specific  enthalpy  for  a  mixture  of  N  species  is  given  by 

h  =  E  cihi 


(53) 
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ht  =  e4  +  RJ  +  (Ai5y)®  (57) 

where  again,  q  is  the  mass  fraction  and  hj  is  the  specific  enthalpy  for  species  i.  The  internal 
energy  at  each  required  temperature  and  pressure  is  calculated  using  Eqs.  (48)-(51)  combined 
with  Eqs.  (54)  and  (55). The  solution  procedure  required  a  previous  call  to  the  equilibrium 
composition  subroutine  to  input  the  partial  pressures  and  mixture  molecular  mass.  The  electronic 
derivatives  appearing  in  Eq.  (51)  are  calculated  using  a  4th-order  four-point  Richardson 
extrapolation  with  AT  =  .OIK.  The  derivatives  were  taken  with  a  frozen  principal  quantum 
number  cutoff.  This  was  done  for  ease  of  calculation  and  to  avoid  the  problem  other  authors14 
experienced  with  large  computational  derivatives  being  obtained  if  the  temperature  and  pressure 
state  corresponded  to  a  partition  function  jump  due  to  the  addition  or  subtraction  of  another 
electronic  level. 


RESULTS  AND  DISCUSSION 

The  chemical  composition  in  terms  of  the  mole  fractions  at  1  atm.  and  Te/T  =  1.0  from  580.2 
K  to  58,020  K  is  shown  in  Fig.  1.  We  see  the  same  rapid  dissociation  of  C2F4  as  observed  by 
Paulino  and  Squires21  due  to  the  inherent  weakness  of  the  carbon  double  bond  in  this  molecule 


as  a  result  of  the  repulsive  electronic  weakening  of  the  bond.  The  results  are  consistent  with 
those  obtained  by  Kovitya.26  Figure  2  shows  a  view  of  the  composition  on  a  logarithmic  scale. 
Figure  3  shows  the  composition  in  the  low  temperature  region  where  molecules  dominate.  In  our 
case  we  see  that  at  about  800  K  C2F4  partially  dissociates  into  C2F2  and  CF4,  recombines 
somewhat  as  temperature  increases,  creating  a  region  where  C2F4,  CF2,  and  CF3  dominate  at 
about  1700  K.  Dissociation  of  C2F4  and  CF2  begins  at  about  2,000  K  into  large  amounts  of  CF3 
which  is  almost  completely  dissociated  by  4,000  K.  The  primary  dissociation  products  are  F  and 
C  which  reach  their  maxima  at  about  5,000  K.  Past  this  temperature  ionization  begins  to  occur 
and  singly  ionized  C  and  F  and  electrons  dominate  the  composition.  Second  ionization  begins 
to  occur  between  22,500  and  25,000  K.  Figures  4  through  9  show  representative  species 
compositions  of  C,  C+,  F,  F+,  e'  and  CF2  at  two  different  pressures,  1  atm  and  .1  atm,  for  the 
case  of  thermal  nonequilibrium.  They  are  given  for  four  heavy  isotherms  of  .05,  .1,1,  and  5 
ev.  Here  we  see  drastically  different  behavior  due  to  the  highly  energetic  electrons  at  their 
elevated  temperature.  Dissociation  of  CF2  happens  rapidly,  other  molecules  exhibited  similar 
behavior.  Single  and  multiple  ionization  of  C  and  F  follows  the  expected  pattern  of  following 
the  electron  temperature,  consistent  with  our  original  formulation.  Figure  8,  the  electron  partial 
pressure,  shows  we  reach  regions  where  dPe/dTe  goes  to  zero.  The  composition  is  constant  after 
this  for  we  had  assumed  that  the  maximum  ionic  charge  state  was  four.  This  suggests  that  for 
a  correct  chemical  model  at  elevated  electron  temperatures  we  may  need  to  extend  the  possible 
ionic  species  to  a  greater  charge  value.  The  degree  if  ionization  is  shown  in  Figs.  10  and  11.  The 
results  show  that  it  tends  to  be  a  stronger  function  of  heavy  particle  temperature  as  pressure 
increases  at  a  given  electron  temperature.  And  reduced  pressures  result  in  greater  degrees  of 
ionization,  as  is  expected.  The  molecular  ionic  and  electronegative  species  were  found  to  exist 
in  only  very  small  amounts  at  all  T/T  values.  Figure  12  shows  the  mixture  enthalpy  and  internal 
energy  for  the  single  temperature  case.  Here  we  see  that  it  is  a  very  strong  function  of 
temperature  throughout  the  molecular  region.  Figures  13  and  14  the  mixture  enthalpy  and 
internal  energy  for  the  two-temperature  case.  We  see  the  most  drastic  changes  with  electron 
temperature  in  the  region  below  10,000  K  in  electron  temperature.  After  this  we  also  see  larger 
differences  with  pressure  for  each  heavy  particle  temperature.  In  future  work  it  is  planned  to 
include  multi-phase  species  in  like  amorphous  carbon.  Molecular  effects  that  are  being  evaluated 
for  inclusion  in  future  models  include  anharmonic  vibrations  and  internal  energy  mode  coupling. 
Also  research  is  being  performed  to  develop  the  complete  reactive  transport  model  for  this 
complex  mixture. 


APPENDIX  1 


For  statistically  independent  subsystems,  the  energy,  entropy  and  other  thermodynamic 
quantities  are  still  additive.  Thus 
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Taking  the  derivatives  of  the  above  we  get 


dS  =Y,  dSn  and  dE  =  E  dEn 


(59) 


Using  the  proper  form  for  the  second  law  for  a  multi-temperature  system5 
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where  5Qn=dEn+5Wn,  and  5Wn=pndV.  Then  we  see  that 
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Now,  for  brevity,  we  invoke  the  relation,  without  proof  that2710 

S  =  kN\n(Q)  +  *p  E 
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And  due  to  the  fact  that  the  n  subsystems  are  statistically  independent  such  that  the  above 
equation,  which  is  derived  for  a  system  in  thermodynamic  equilibrium  also  holds  for  all  separate 
subsystems  then  we  write5 

Sn  =kNhx(Qn)  En  (63) 

so  that 

(dS)NJy  =  kfidE  (64) 


therefore 
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For  our  two-temperature  system 
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thus 
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Now  we  define  an  "economy"  between  the  heavy  particle  and  electron  energies,  q,  which  defines 
the  fractional  contribution  of  each  to  the  total  energy  in  any  state.  Lets  express  it  as 


Ee  =  r)E 


(69) 


and 
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Then 
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and 


dEH  =  (1  -  r\)dE 


(72) 


which  upon  substitution  into  Eq.  (68),  and  the  subsequent  cancellation  of  dE  yields  that 
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with  the  limits  that  if 
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Figure  1.  Teflon  composition  at  1  atm.  for  isothermal  case. 
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Figure  3.  Isothermal  Teflon  composition  for  T  <  5,000  K. 
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Figure  4.  Partial  pressure  of  C  vs.  Te. 
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Figure  5.  Partial  pressure  of  C+  vs.  Te. 
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Figure  6.  Partial  pressure  of  F  vs.  Te. 
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Figure  7.  Partial  pressure  of  F+  vs.  Te. 
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Figure  8.  Partial  pressure  of  electrons  vs.  Te. 
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Figure  9.  Partial  pressure  of  CF2  vs.  Te. 
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Figure  10.  Degree  of  Ionization  vs.  Te 
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Figure  11.  Degree  of  Ionization  vs.  Te,  log  scale. 


Figure  12.  Single  temperature  mixture  enthalpy  and  internal 
energy  vs.  T  at  1  atm. 
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Figure  13.  Mixture  enthalpy  vs.  Te. 
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Figure  14.  Mixture  internal  energy  vs.  Te. 
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ABSTRACT 

Electron-beam  diodes  driven  by  fast-rising,  high-voltage  pulses  often  operate 
with  cold  cathodes  for  which  the  presence  of  a  plasma  adjacent  to  the  cathode  surface  is 
essential  to  obtain  adequate  electron  emission.  A  consequence  of  such  surface  plasma, 
however,  is  closure  of  the  interelectrode  gap  by  plasma  motion.  The  diode  impedance 
decreases  with  time,  adversely  affecting  the  efficiency  of  coupling  to  the  power  source. 

Plasma  closure  of  the  diode  gap  also  limits  the  length  of  the  electron  beam  pulse,  and 
the  ability  to  operate  the  diode  repetitively  at  high  frequency.  Resistive  heating  of  the 
plasma  competes  with  work  performed  in  expanding  the  plasma  and  heat  transfer  to  the 
cold-cathode  boundary.  The  resulting  closure  speed  is  calculated,  using  an  MHD  code, 
and  found  to  agree  well  with  results  of  experiments  using  organic-cloth  cathodes  at  35 
kV.  Computed  plasma  speeds  are  typically  8-12  km/s,  and  are  relatively  insensitive  to 
the  applied  voltage.  Gap  closure  due  to  the  plasma  motion  calculated  numerically 
corresponds  to  estimates  based  on  impedance  collapse  in  the  experiments. 

INTRODUCTION 

Closure  of  the  interelectrode  space  in  electron-beam  diodes  operating  with  cold-cathodes  is  a 
well-known  phenomena  limiting  the  performance  of  such  devices.  Plasma  adjacent  to  the  cathode 
surface  is  needed  to  obtain  adequate  electron  emission  for  space-charge  limited  flow  in  the  diode.  Such 
plasma,  however,  can  cross  the  diode  gap,  resulting  in  collapse  of  the  diode  impedance  during  the  high- 
voltage  pulse.  This  behavior  adversely  affects  the  efficiency  of  coupling  between  the  electron  beam  and 
the  power  source,  limits  the  duration  of  the  electron-beam  pulse,  and  may  preclude  operation  at  desired 
repetition  rates.  It  is  useful,  therefore,  to  understand  factors  that  determine  plasma  closure  in  electron- 
beam  diodes  in  the  context  of  theoretical  modeling  that  can  then  be  employed  to  examine  directions  for 
improvement.  The  first  step  in  such  modeling  consists  of  simulating  existing  experiments  using 
techniques  that  can  be  extended  to  more  complex  possibilities. 

The  simplest  representation  of  plasma  closure  in  a  diode  comprises  a  thin  layer  of  uniform 
plasma  between  a  solid  cathode  surface  and  a  region  representing  the  vacuum  in  the  diode  gap.  The 
physical  model  consists  of  expansion  of  the  plasma  across  the  gap,  with  resistive  heating  compensating 
for  the  loss  of  internal  energy  to  work  and  heat  transfer.  Even  this  elementary  model  is  too  complex  for 
accurate  analytical  treatment,  so  numerical  techniques  are  necessary.  Estimates  of  particle  density  and 
temperature  for  the  surface  plasma  suggest  that  a  continuum  approach  would  be  appropriate  for  the 
motion  of  the  main  portion  of  the  plasma.  The  MACH2  2-1/2  dimensional  MHD  code  is  therefore  a 
reasonable  choice  of  calculational  tool.  Such  a  code  is,  of  course,  clearly  inadequate  to  model  the 
behavior  of  an  electrically  non-neutral,  non-continuum  flow. 

The  essential  feature,  however,  that  must  be  calculated  in  order  to  compute  the  dynamics  and 
thermodynamics  of  the  surface  plasma  is  the  current  density.  For  MACH2,  this  merely  requires  that  the 
resistivity  of  the  “vacuum’  region  has  a  value  corresponding  to  space-charge  limited  flow  for  the 
instantaneous  values  of  applied  voltage  and  vacuum  gap.  For  most  situations,  this  resistivity  is 
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sufficiently  high  to  permit  magnetio-flux  diffusion  rapidly  through  the  “vacuum’  region,  which  results  in 
nearly  uniform  current  density  normal  to  the  surface  of  the  plasma  layer.  Thermal  expansion  of  the 
plasma  decreases  the  effective  diode  gap,  but  the  resistivity  of  the  vacuum  region  may  be  continually 
adjusted  to  maintain  the  proper  value  of  space-charge  limited  current  density  during  the  voltage  pulse. 
While  this  technique  cannot  simulate  the  actual  behavior  of  the  particle  flow  in  the  vacuum,  it  is 
adequate  for  providing  the  correct  power  to  the  plasma  in  order  to  compute  the  motion  of  the  main 
portion  of  the  plasma  mass.  Detailed  consideration  may  then  be  given  to  penetration  of  electric  field  into 
the  plasma  in  defining  the  position  of  the  electron  emitting  surface.  This  is  accomplished  in  terms  of  the 
local  values  of  the  Debye  length  at  the  vacuum  edge  of  the  plasma,  aided  by  the  condition  of  electric 
field  equal  to  zero  for  emission  into  space-charge  limited  flow. 

USE  OF  THE  MACH2  CODE 

The  MACH2  code  [1]  allows  numerical  solution  of  the  equations  for  2-1/2  dimensional,  unsteady, 
MHD  flow.  It  uses  an  arbitrary  Lagrangian-Eulerian  (ALE)  formulation  in  combination  with  a  multiple 
block  format  that  permits  application  to  complex  geometries.  In  addition  to  classical  transport  properties 
and  ideal  gas  formulas,  MACH2  operates  with  the  SESAME  tables  for  equation-of-state  and  transport 
properties  of  a  variety  of  materials.  Anomalous  resistivity  coefficients  based  on  plasma/beam 
microinstabilities  are  also  available,  as  are  options  for  optically-thin  and  both  equilibrium  and  non¬ 
equilibrium  diffusion  models  for  radiative  transport.  Several  opportunities  exist  as  well  for  coupling 
plasma  dynamics  to  external  electrical  sources,  including  prescribed  voltage  or  current  waveforms  and 
self-consistent  solution  of  lumped-element  circuitry. 

MACH2  is  nevertheless  a  continuum  formulation,  based  on  electrical  quasi-neutrality  and  local 
(two-  or  three-  temperature)  thermodynamic  equilibrium.  Its  application  to  problems  in  an  electron-beam 
diode  is  restricted  to  regions  where  the  plasma  density  is  high  enough  to  satisfy  these  physical 
constraints.  Plasma  created  at  the  cathode  surface  by  explosion  of  asperities  at  the  cathode  surface,  or 
by  injection  from  a  continuum  plasma  source,  may  satisfy  the  continuum,  quasi-neutral,  LTE  conditions 
adequately  for  the  purpose  of  calculating  plasma  motion.  The  very  low  density,  non-neutral  region 
between  the  emitting  surface  and  the  anode  certainly  does  not.  MACH2,  however,  permits  the  electrical 
resistivity  of  a  region  or  material  (defined,  for  example,  by  its  low  mass-density)  to  be  assigned  quite 
flexibly.  Thus,  very  low-density  material  in  the  vacuum  gap  can  be  given,  at  any  instant  in  the 
calculation,  a  value  of  resistivity  that  will  provide  the  correct  total  impedance  corresponding  to  space- 
charge  limited  flow.  For  the  applied  voltage,  therefore,  the  correct  total  current  will  be  delivered  through 
the  expanding  surface-plasma.  Implicit  in  this  approach  is  the  approximation  that  the  basically  one¬ 
dimensional  formulation  of  space-charge  limited  flow  will  still  be  adequate  to  describe  operation  of  the 
vacum  region.  To  correct  the  value  of  the  current  from  the  Child-Langmuir  formula,  in  order  to  obtain  the 
actual  experimental  current  before  significant  impedance  collapse  has  occurred  due  to  plasma  motion,  a 
value  of  perveance  is  calculated.  This  correction  factor  is  then  held  constant  during  diode  closure. 

MACH2  has  not  been  used  to  calculate  the  actual  explosion  of  surface  asperities  to  create  the 
initial  plasma.  Instead,  the  conditions  of  this  initial  plasma  (mass  density,  electron  and  heavy-particle 
temperatures,  and  initial  thickness)  are  assumed  to  be  uniform,  and  treated  as  parameters  in  analyzing 
the  subsequent  plasma  expansion.  For  properties,  such  as  the  closure  speed,  non-dimensionalization 
suggests  that  the  the  actual  value  of  initial  mass-density  and  initial  plasma  thickness  are  not  crucial. 
Similarly,  low  estimates  of  the  initial  temperature  might  be  self-correcting  due  to  enhanced  resistive 
heating.  (Over-estimates  of  the  temperature,  however,  will  simply  cause  the  plasma  to  expand  too 
quickly  in  comparison  to  experimental  data.) 

SPECIFICATION  OF  THE  ELECTRON-EMISSION  SURFACE 

Before  proceeding  to  solve  the  dynamics  of  plasma  closure,  and  to  examine  the  sensitivities  of 
the  model  to  the  choice  of  initial  plasma  conditions,  it  is  necessary  to  develop  a  specification  for  the 
position  of  the  electron-emitting  surface.  For  the  simple  case  of  emission  from  a  solid  surface,  the  diode 
gap  is  merely  the  geometric  distance  between  the  electrodes,  to  the  level  of  atomic  dimensions. 
Emission  from  a  plasma  surface  is  less  well-defined  and  can  involve  uncertainties  in  effective  diode-gap 


on  the  order  of  the  gradient  scale-length  for  the  plasma  density,  the  mean  free  path  for  electron-electron 
collisions,  the  (local)  Debye  length  in  the  plasma,  or  some  other  condition.  For  the  present  calculations, 
we  have  chosen  to  define  the  electron-emitting  surface  of  the  plasma  based  on  the  condition  that  the 
electron  velocity-distribution  is  significantly  distorted  from  a  Maxwellian  distribution  due  to  the  electric 
field  near  the  emitting  surface.  Algebraically,  this  condition  is  written  as: 


E(x)  Xd  >  kTe 


(D 


where  Xq  is  the  Debye  length  based  on  the  local  plasma  temperature  Te  and  the  local  electron 
density,  which  varies  sharply  from  the  plasma  layer  to  the  vacuum  region.  The  local  electric  field,  E(x), 
would  be  zero  exactly  at  the  emitting  surface,  as  a  condition  of  space-charge  limited  flow,  but  has  a  finite 
value  at  a  distance  of  a  Debye  length  into  diode  gap,  i.e.,  x  =  Xo  .  The  condition  of  Eqn.  1  is  thus 
satisfied,  for  any  plasma-electron  temperature,  when  the  electron  density  falls  below  some  value.  While 
the  condition  could  equally  well  have  been  written  with  a  factor  (of  order  unity)  multiplying  the  electron 
temperature,  the  resulting  variations  in  diode  gap  would  usually  be  less  than  the  resolution  of  the 
calculational  grid  employed  in  the  present  work.  Furthermore,  within  the  approximation  of  using  the  one¬ 
dimensional,  Child-Langmuir  formula  for  space-charged  limited  flow,  we  monitor  the  diode  gap  according 
to  Eqn.  1  only  at  one  position  along  the  surface  (e.g.,  half-radius).  A  more  laborious  approach  involving 
local  calculation  of  the  diode  gap  is  deferred  until  warranted  by  substantially  two-dimensional  motion  of 
the  plasma  surface  and  until  an  appropriate  substitute  for  the  Child-Langmuir  formula  can  be  developed, 
(which  probably  will  require  a  particle-in-cell  calculation). 

CALCULATIONAL  RESULTS 

The  first  application  of  the  present  model  has  been  to  a  set  of  experiments  [2]  performed  with 
carbon-cloth  cathodes  at  relatively  modest  applied  voltages  (~  35  kV).  These  experiments  were  chosen 
for  simulation  because  the  clean  and  simple  voltage  waveform  (linear  rise  and  fall  surrounding  a 
relatively  long  constant  value)  allowed  less  uncertainty  in  comparing  theoretical  and  experimental 
behavior.  The  experimental  arrangement  comprises  a  cylindrical  cathode,  2  cm  in  radius,  with  a  circular 
endface  covered  by  the  emitting  carbon-cloth,  separated  by  a  gap  of  0.5  cm  from  a  screen  anode;  (the 
experiments  were  directed  toward  virtual-cathode  operation  for  microwave  production.)  Figure  1  displays 
the  experimental  and  theoretical  waveforms  for  voltage  and  current.  As  previously  noted,  a  perveance 
factor  is  included  in  the  calculation  to  match  space-charge  limited  flow  to  the  experimental  current  value 
before  significant  diode  closure  has  occurred  (t  <  50  ns).  This  factor  is  held  constant  during  the 
subsequent  plasma  motion.  The  computed  current  waveform  agrees  rather  well  with  the  experimental 
curve  scanned  from  the  published  oscillogram.  Experimental  behavior  from  210  to  300  ns  is  difficult  to 
discern  due  to  relatively  high  amplitude,  high-frequency  oscillations  on  the  trace;  the  calculations  assume 
an  inductance  of  100  nh  in  series  with  the  diode,  which  may  be  somewhat  too  high. 

Figure  2  provides  the  calculated  voltage  and  current  waveforms  along  with  the  diode  gap.  In 
these  calculations,  which  serve  as  a  baseline  for  parameter  variations,  the  initial  electron  and  ion 
temperatures  in  the  plasma  layer  are  1.0  eV  and  0.1  eV,  respectively,  and  the  initial  plasma  (mass) 
density  is  10'3  kg/m3.  (These  values  are  chosen  arbitrarily,  but  should  be  relevant  to  low  temperature 
plasmas.)  The  material  of  the  plasma  is  carbon  phenolic,  available  in  the  SESAME  tables,  and  was 
chosen  in  order  to  account  for  the  complicated  molecular  composition  of  the  carbon-cloth  cathode,  which 
consists  of  silk  and  rayon  (vs  pure  carbon).  The  effective  molecular  mass  of  the  phenolic  is  9.01  amu. 
The  initial  thickness  of  the  surface  plasma  is  0.5  mm,  but  rapidly  increases  to  about  1  mm,  due  to  the 
application.#  Eqn.  1  to  the  density  distribution  computed  at  early  times. 

The  speed  of  plasma  closure  is  not  constant,  but  corresponds  to  an  average  speed  of  8  km/s. 
This  value  is  consistent  with  the  estimates  from  the  experimentally-observed  decrease  of  diode 
impedance  with  time  (assuming  Child-Langmuir  behavior).  Such  consistency,  of  course,  is  directly 
associated  with  the  good  agreement  between  experiment  and  theory  for  the  diode  current.  In  addition  to 
the  overall  behavior  of  diode  impedance  -collapse  due  to  plasma  motion,  the  MACH2  simulations 


provide  detailed  descriptions  of  plasma  properties  and  distributions  within  the  diode.  Examples  of  this 
information  are  given  in  Fig.  3,  which  displays  the  shape  of  the  interface  between  plasma  and  vacuum, 
and  contours  of  mass  density  and  electron  temperature. 

SENSITIVITIES  TO  PARAMETER  CHOICES 

The  choices  of  initial  values  for  the  plasma  layer  may  be  informed  by  experience,  but  are 
certainly  arbitrary  in  the  present  work.  It  is  necessary,  therefore,  to  explore  the  effects  of  other  values  on 
the  basic  behavior  of  diode  closure.  Three  important  parameters  are  the  initial  values  of  electron 
temperature,  ion  temperature  and  mass  density.  In  the  following  comparisons,  the  voltage  waveform  will 
be  maintained,  and  only  one  of  these  parameters  will  be  varied  at  a  time;  the  other  parameters  will  be 
held  at  their  values  for  the  baseline  case  previously  discussed. 

Figure  4  displays  the  current  and  diode  gap  histories  for  initial  values  of  electron  temperature 
that  are  factors  of  two  lower  and  higher  than  the  baseline  case,  0.5  eV  and  2.0  eV,  respectively.  The 
lower  value  results  in  little  change  in  diode  gap  and  corresponding  small  change  in  the  diode  impedance 
with  time;  the  current  follows  the  voltage  waveform.  On  the  other  hand,  the  higher  value  of  initial  electron 
temperature  results  in  substantial  faster  diode  closure  and  a  much  higher  value  of  diode  current  than  is 
observed  experimentally.  Comparison  of  the  diode  gap  histories  with  that  of  the  baseline  case  indicates 
the  importance  of  the  early  development  of  closure  speed  on  subsequent  decrease  of  the  diode  gap  and 
rapid  increase  of  diode  current  (which  scales  inversely  with  the  square  of  the  gap).  Higher  current  density 
helps  to  maintain  the  plasma  internal-energy  as  the  plasma  expands  across  the  gap.The  early 
development  of  closure  speed  in  the  higher  temperature  case  is  simply  a  result  of  higher  (electron) 
pressure. 

In  Fig.  5,  the  initial  value  of  electron  temperature  is  1.0  eV,  but  the  initial  values  of  ion 
temperature  have  been  increased  by  factors  of  five  and  ten  over  the  baseline  case,  0.5  eV  and  1 .0  eV, 
respectively.  Higher  initial  temperatures  for  the  ions  also  means  higher  plasma  pressure  (at  fixed  mass 
density),  so  closure  speeds  develop  faster  than  in  the  baseline  case.  There  is  not  as  much  sensitivity, 
however,  to  variation  in  ion  temperature  (vs  initial  electron  temperature)  because  the  ions  in  the  baseline 
case  rapidly  warm  up  to  above  0.5  eV  due  to  heat  transfer  from  the  electrons.  The  effect  of  changing  the 
initial  value  of  mass  density  is  much  less  pronounced  because,  at  fixed  initial  values  of  electron  and  ion 
temperatures,  the  plasma  pressure  scales  approximately  with  the  mass  density;  (variations  are  possible 
due  to  changes  in  ionization  with  density).  Thus,  the  basic  expansion  speed  of  the  plasma  remains 
approximately  the  same  as  for  the  baseline  case.  The  results  of  this  are  observed  in  Fig.  6  for  which 
factors  of  ten  higher  and  lower  values  of  initial  mass  density  are  used. 

To  the  extent  that  the  simulations  have  captured  the  actual  behavior  of  diode  closure,  the 
agreement  of  the  baseline  case  with  the  experimental  data  suggests  that  the  properties  of  the  initial 
plasma  are  reasonably  given  by  electron  and  ion  temperatures  of  1  eV  and  0.1  eV,  respectively.  The 
initial  mass  density,  however,  could  be  higher  or  lower  by  almost  a  factor  of  ten  without  disturbing  this 
agreement  very  much. 


EFFECTS  OF  APPLIED  VOLTAGE 

The  present  model  can  be  employed  to  examine  the  effects  of  variation  of  circuit  operation  on 
diode  closure.  For  example,  if  the  same  temporal  waveform  is  provided  to  the  diode  with  different 
amplitudes,  the  change  in  diode  closure  (and  thus  diode  impedance  history)  can  be  predicted.  Figure  7 
displays  the  calculated  diode  gap  histories  as  the  amplitude  of  the  applied  voltage  is  increased  by  factors 
of  1 .5  -  3  over  the  baseline  case.  The  average  closure  speed  only  varies  by  a  factor  of  50  %  ,  which 
agrees  with  experience  in  a  variety  of  cold-cathode  diodes.  The  peak  speed,  however,  shows  greater 
sensitivity  to  applied  voltage.  Examination  of  the  two-dimensional  distributions  of  plasma  in  the  diode 
indicate  that  greater  nonuniformity,  especially  faster  closure  near  the  centerline  of  the  diode  may  be 
responsible  for  these  higher  speeds  elsewhere  in  the  diode;  (recall  that  the  diode  gap  is  only  monitored 
at  mid-radius  in  these  calculations).  The  basic  result  of  the  relative  insensitivity  of  diode  impedance- 


collapse  to  applied  voltage  is,  perhaps,  the  only  claim  that  may  be  made  until  the  present  calculations 
are  extended  to  more  two-dimensional  treatment  of  current  flow  in  the  vacuum. 

CONCLUDING  REMARKS 

It  appears  that  MHD  codes,  such  as  MACH2,  can  be  successfully  employed  to  study  plasma 
phenomena  in  high  current-density  diodes.  The  basis  for  this  success  at  present  is  the  application  to 
diode  problems  for  which  one-dimensional  approximations  to  the  current  flow  in  the  vacuum  are 
adequate.  Also,  by  invoking  the  condition  of  space-charge  limited  flow,  (E(0)  ~  0),  the  specification  of  the 
electron-emitting  surface  is  facilitated.  For  two-  or  three-dimensional  problems,  recourse  to  particle-in- 
cell  (PIC)  techniques  to  calculate  the  local  current  density  through  the  plasma  layer  will  probably  be 
necessary.  For  some  problems,  the  development  of  the  current  flow  and  charge  distribution  in  the 
vacuum  may  occur  on  a  timescale  much  shorter  than  the  hydrodynamic  expansion  times  for  the  surface 
plasma.  This  situation  may  then  allow  calculation  of  diode  behavior  with  surface  plasmas  by  a  rapid 
iteration  between  PIC  codes,  run  for  several  nanoseconds,  and  hydrocodes  running  over  a  microsecond, 
without  encountering  the  difficulties  of  including  plasma  motion  within  the  PIC  code  itself. 
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Fig.  3.  MACH2 
provides  2D 
distributions  of 
plasma  properties 
including  plasma- 
vacuum  interface. 


electron 

temperature  (eV)  —  3.2082E-01 


+-  1 . 4413E+00 


0  0  05  O.t  0.15  02  0  2S  02  0 

« _ tjmefusl _  , 


0.1  015  0.2  025  0.: 

time  (ms)  _ 


0  0.05  O.t  0.15  0.2  025  0.3  0  005  0.1  0.15  0.2  025  0.: 

time  (ps)  -j-e  ^  2  time  (ps) 

Fig.  4.  Other  initial  electron  temperatures  (0.5 
and  2.0  eV  vs.  1 .0  eV)  predict  plasma  closure 
that  is  inconsistent  with  experiment. 
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Fig.  6.  Other  initial  mass  densities  at  the  surface 
(1  O'2  and  1 0'4  kg/m3  vs.  10'3  kg/m3)  predict  closure 
that  is  less  consistent  with  experiment. 
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Fig.  5.  Other  initial  ion  temperatures  (0.5  and 
1 .0  eV  vs.  0.1  eV)  also  predict  closure  that  is 
inconsistent  with  experiment. 
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Fig.  7.  The  time  dependent  gap  is  computed 
for  different  applied  voltages  relative  to  the  35 
kV  baseline  case.  The  average  speed  is 
between  8  and  12  km/s. 


